Problem in verifying zero blocking property of a zero of a multiple input multiple output system

The transmission zeros of a MIMO system are obtained by solving

\begin{align}
\dot{x}^* = A\,x^* + B\,u^*, \\
y = C\,x^* + D\,u^*,
\end{align}

for $x^*$ and $u^*$ such that $y = 0$. By using the Laplace variable $s$ for the time derivative that system of equation can also be written as

$$
\begin{bmatrix}
0 \\ 0
\end{bmatrix} =
\begin{bmatrix}
A - s\,I & B \\ C & D
\end{bmatrix}
\begin{bmatrix}
x^* \\ u^*
\end{bmatrix},
$$

which can be seen as the generalized eigenvalue problem $\mathcal{A}\,v = \lambda\,\mathcal{B}\,v$, with $\lambda$ equivalent to $s$,

$$
\mathcal{A} =
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}, \quad
\mathcal{B} =
\begin{bmatrix}
I & 0 \\ 0 & 0
\end{bmatrix}, \quad
v =
\begin{bmatrix}
x^* \\ u^*
\end{bmatrix}.
$$

I believe that tzero in matlab also solves such generalized eigenvalue problem to obtain the transmission zeros (however by using eig yourself the associated eigenvectors are more accurate than the nullspace calculations). The associated solution in time for both $x(t)$ and $u(t)$ can be obtained using

$$
\begin{bmatrix}
x(t) \\ u(t)
\end{bmatrix} = e^{\lambda\,t}\,
\begin{bmatrix}
x^* \\ u^*
\end{bmatrix},
$$

so you can't just pick any value for $\lambda$ when generating such solution. You could now simulate the system using $x(0)=x^*$ as initial condition and the expression for $u(t) = e^{\lambda\,t}u^*$ with lsim. However, I believe that slim uses zero-order-hold for the input so the actual applied input will deviate slightly causing the output $y$ to deviate away from zero using that approach. Instead, if you look at the analytical solution, both the output and the error in state dynamics should remain practically to zero (for non-positive $\lambda$):

  n = length(A);
  m = length(D);
  [V,L] = eig([A B; C D], diag([ones(1,n) zeros(1,m)]));
  ind = 1;    % Index of associated eigenvalue (avoid infinite eigenvalues)
  lambda = L(ind,ind);
  x0 = V(1:n, ind);
  u0 = V(1+n:n+m, ind);

  N = 1e3;
  t = linspace(0, 5, N);
  x = @(t) x0 * exp(lambda* t);
  u = @(t) u0 * exp(lambda* t);
  y = @(t) C * x(t) + D * u(t);
  e = @(t) A * x(t) + B * u(t) - lambda* x(t);

  Y = zeros(n, N);
  E = zeros(m, N);
  for k = 1 : N
      Y(:,k) = y(t(k));
      E(:,k) = e(t(k));
  end
  figure, plot(t, E)
  figure, plot(t, Y)


Why call it discrete time signals as difference equations?

A difference equation is a recurrence relation where the current element $y[n]$ of a sequence is related to its past values $y[n-k]$, $k>0$ (if we interpret the index $n$ as time index, which is not necessary). The equation you gave in your question is a special case, because there are no past terms $y[n-k]$ involved, so you could call it a zeroth order difference equation, but you're right that from a common sense point of view it is no difference equation at all.

However, what counts here is the fact that no past (or future) elements $y[n-k]$ are involved, and not that "nothing is getting subtracted". E.g., the following equation is a standard example of a difference equation:

$$y[n]=ay[n-1]+bx[n]\tag{1}$$

The subtraction takes place in the index $n$ ($y[n-1]$). But obviously you can substitute $n+1$ for $n$ yielding

$$y[n+1]=ay[n]+bx[n+1]\tag{2}$$

which is completely equivalent to $(1)$ with no subtraction taking place, but it's still a difference equation, because the difference between the indices remains unchanged.


Not really, as the transform is real. However, one could interpret the sign as a poor man's phase, being "quantized" or restricted to values $0$ or $\pi$. In other words, $1 = 1.e^{0.\imath}$ and $-1 = 1.e^{\pi.\imath}$.

[EDIT] There are some instances where people use the sign of DCT or (real) wavelet coefficients, for subpixel image registration or matching, as alluded to in DCT Sign Only Correlation and Its Application to Image Registration:

  In the DCT case, all DCT coefficients are real values, as a result,
  the phase information corresponds to the sign of DCT coefficients.


downsample signal to a specific number of samples in Python

I am doing something similar to your application with 1D CNN. I think scipy.resample_poly is the most versatile function since it allows both upsampling, downsampling, or a combination of both.

Also your comment about "values beyond the boundary of the signal to be zero" can be solved by using the option "line" in padtype, as shown in the function documentation, in order to eliminate spurious effects at the head/tail of the signal.

In general if the integer number of resampling factor is still an issue to you (but it shouldn't as you can play with the length of the signal) you can play with a combination of moving average followed by interpolation on the grid.


I have some questions about background subtraction with Gaussian Mixture Model (GMM)

Yes, the dimension of $X$ should be 3, one for each of the red, green, and blue channels. You can see this from their Figure 2, where they show plots of the distributions of red and green channels only (to motivate why they want to use a Gaussian Mixture Model in the first place --- the distribute of pixel values is multi-modal).
Yes, $I$ is the identity matrix.
Random variables can have the same variances and still be independent. The variance is part of the model. In this case what the model is saying is: the red, green, and blue pixel values vary without regard to each other (they are independent) but the amount that they vary is about the same as each other (they have the same variances).
Whether only 1 distribution is chosen when $K=3$ will depend on the threshold $T$, surely?  Yes, it is possible that only one distribution will be chosen.


We could say MUSIC is optimal algorithm for spectrum estimation if the number of spectral peaks is a priory known and is less than estimator order, $M$. In MUSIC we need to estimate $R_{xx}$ autocovariance matrix of input random process first. The matrix dimension is limited to estimator's order, $M$. Second step is to perform SVD: $R_{xx} = U\cdot\Sigma\cdot V^{H}$ and $\Sigma$ is diagonal matrix with mimension also $M$. This matrix contains singular values of $R_{xx}$, the largest $K$ values correspond to signal's subspace while other - to noise one. For obvious reason $K < M$, so the limitation of the number of spectral peaks (in the other words - harmonic signal sources) becomes clear.

When you use MUSIC in practice to estimate PSD you need either to know actual number of not noise signals in the mixture or to estimate this number after SVD by some kind of threshold or ML search. If your a priory information or estimation isn't valid, it can strongly violate algorithm performance.

What you're searching is number of sources estimation or source enumeration. You can look through this article for example. In the case of MUSIC source enumeration problem can be viwed as signal's subspace estimation problem. You can find simple approach for the problem's solving in this article. It is called MDL, I checked it once while modeling the described SNR estimation algorithm in MATLAB. It seems to work well for pretty low SNR values. But remember you can't overcome limitation of $K < M$ for subspace based methods as I've described above.

Hope it helps.


How do I implement cross-correlation to prove two audio files are similar?

Cross-correlation and convolution are closely related.  In short, to do convolution with FFTs, you

zero-pad the input signals (add zeros to the end so that at least half of the wave is "blank")
take the FFT of both signals
multiply the results together (element-wise multiplication)
do the inverse FFT

conv(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros))

You need to do the zero-padding because the FFT method is actually circular cross-correlation, meaning the signal wraps around at the ends.  So you add enough zeros to get rid of the overlap, to simulate a signal that is zero out to infinity.

To get cross-correlation instead of convolution, you either need to time-reverse one of the signals before doing the FFT, or take the complex conjugate of one of the signals after the FFT:

corr(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros[reversed]))
corr(a, b) = ifft(fft(a_and_zeros) * conj(fft(b_and_zeros)))

whichever is easier with your hardware/software.  For autocorrelation (cross-correlation of a signal with itself), it's better to do the complex conjugate, because then you only need to calculate the FFT once.

If the signals are real, you can use real FFTs (RFFT/IRFFT) and save half your computation time by only calculating half of the spectrum.

Also you can save computation time by padding to a larger size that the FFT is optimized for (such as a 5-smooth number for FFTPACK, a ~13-smooth number for FFTW, or a power of 2 for a simple hardware implementation).

Here's an example in Python of FFT correlation compared with brute-force correlation:  https://stackoverflow.com/a/1768140/125507

This will give you the cross-correlation function, which is a measure of similarity vs offset.  To get the offset at which the waves are "lined up" with each other, there will be a peak in the correlation function:

The x value of the peak is the offset, which could be negative or positive.

I've only seen this used to find the offset between two waves.  You can get a more precise estimate of the offset (better than the resolution of your samples) by using parabolic/quadratic interpolation on the peak.

To get a similarity value between -1 and 1 (a negative value indicating one of the signals decreases as the other increases) you'd need to scale the amplitude according to the length of the inputs, length of the FFT, your particular FFT implementation's scaling, etc.  The autocorrelation of a wave with itself will give you the value of the maximum possible match.

Note that this will only work on waves that have the same shape.  If they've been sampled on different hardware or have some noise added, but otherwise still have the same shape, this comparison will work, but if the wave shape has been changed by filtering or phase shifts, they may sound the same, but won't correlate as well.


Pixel time series trend: What are the causes?

You are very likely to observe trends in either your camera or your light source.

Sources from the camera can include:

thermal effects inside the sensor (basically, the more you use it, the more current will traverse it in various parts and it will get warmer)
drift in the ADC of your sensor (probably linked to the aforementioned).

Also, light sources can also drift with time: thermal effects will influence light bulbs.

You can try to mitigate these effects by:

changing your light source (I would assume LED lamps are more stable, but this needs to be checked)
waiting long enough before starting the measurements in order to avoid transitory phenomenon
carefully playing with the settings of your camera. Typically, I would fix manually all the exposure-related settings (exposure time, gain, framerate...) so that the autocontrols do not make vary the exposure values when an object enters the scene. However, this will prevent the camera from tracking natural scene trends.
simply fitting a linear regression to your signal in order to remove the trend while processing.


Difference between 'conventional' and 'adaptive' beamformers?

A beamformer is basically a spatial filter.  It can be passive, just like a temporal filter.

Instead of samples separated by time, they are separated by space.  A passive temporal filter can be a bandpass that is "aimed" or "steered" at a particular frequency.  For passive spatial filters (i.e. beamformers), the filter can be steered towards a particular angle of arrival, instead of temporal frequency.

Adaptive filters/beamformers can be incredible because they can "steer" a null towards a frequency or angle where there is an interfering signal.

Just like a temporal adaptive filter, a spatial adaptive filter (i.e. beamformer) is constantly adjusting the filter weights/coefficients to optimize some criteria that usually involves "nulling out" or "rejecting" an interferer.

Here's a diagram of a CBF that is called a "K-Omega" Beamformer.

Here's a diagram to reinforce the idea that a passive beamformer is possible.

I realize that these are a bit random, but hopefully you can follow the logic of the above.
I'll see if I can find a better diagram that more clearly shows what's happening.  To be clear, the windows referred to above are Hamming/Hanning type windows, and this process is basically returning a 2-d matrix where temporal frequency is on the x-axis and a special spatial variable is on the y-axis.  This special variables makes the math easier, and it's a one step process to convert the special spatial variable and temporal frequency to an angle of arrival.

The green lines below are lines of constant angles.

Keep in mind that all of this is for CBF (conventional beamformers) and the above covers all spatial and temporal frequency (within the limits of Nyquist.
Some common ABF techniques are:

MPDR - Minimum Power Distortionless Rejection BF (Described by Van Trees)
DMR - Dominant Mode Rejection BF (Abraham and Owsley)
R-DMR - Robust Dominant Mode Rejection BF (Cox and Pitre)
EBAE - DMR BF with Eigenvector Beam Association and Excision (Kogon)

Instead of using a spatial FFT as the spatial filtering operation, these techniques typically involve constructing an updated spatial autocorrelation matrix of the incoming signal, and then using that matrix or eigenvectors of that matrix to adaptively influence the spatial filter.

Update for @Mohammad:
Below is the list that of beamforming texts I received from my professor:

Van Veen, B.D.; Buckley, K.M.; , "Beamforming: a versatile approach to spatial filtering," ASSP Magazine, IEEE , vol.5, no.2, pp.4-24, April 1988
Efficient digital beamforming in the frequency domain
Brian Maranda, J. Acoust. Soc. Am. 86, 1813 (1989)
Array Signal Processing by Don Johnson and Dan Dudgeon, Prentice Hall, 1993
Optimum Array Processing by Harry L. Van Trees, Wiley, 2002.

Also, after searching around, this looks very interesting.  It seems to be more interested in the practical/implementation side than the theory.  I don't have it but I'll probably buy a copy:

Practical Array Processing by Mark Sullivan, McGraw Hill, 2008 ( http://www.practicalarrayprocessing.com/index.html )


Image datasets used for salient objects?

Datasets explicitly mentioning saliency

MSRA
SOD based on BSD Berkeley segmentation dataset

Apparently more general datasets (not researched in detail)

DUT-Omron Image dataset
THUS
ECSD
PASCAL
...

More challenging situations with more noise

Detection on magnetism such as MR: Medical image dataset and Extract voice from magnetic tape


Probably no, but... maybe. The classical sampling theorem is often misstated: if a continuous signal is band-limited, it can be recovered exactly (in theory) if the sampling frequency is twice its maximum frequency.
However, there are cases where signals can be recovered above the Nyquist limit, provided that we have more information, because ambiguities may happen  in the recovery:

finite-time signal cannot be strictly band-limited, but processing can afford a small amount of distorsion,
some band-limited signals can be recovered with about twice their bandwidth (and not just twice the max frequency).
there are more generic sampling theorems (like: do you sample the derivative of the signal at the same time), but more involved as well.

Generally, the above possibilities are not implemented in traditional hardware or software.


The idea behind nonuniform quantization (for example an a, or mu law for telephony speech) is to make better use of the input signal amplitude distribution in lossy source coding applications..
Over a period of time, it's emprically observed that, most of the speech occurs at low amplitude levels, and occasionally at high amplitude (full dynamic range) levels.
When this is the case, a uniform quantizer would be wasting its bits (its quantization levels) into many of those statistically less used louder signal levels.
In order to make efficient use of the quantizer (number of levels and bits), and to minimize the average quantization error (the distortion), it will be wiser to allocate more bits (finer levels) to the (often used) lower signal levels, while assigning less bits, and (coarse levels) to the (less used) louder signal levels.
The technical term for such a modification to a uniform quantizer is called a pdf optimized Max-Lloyd quantizer. The compander is a simplified example by assuming a fixed typical pdf of input speech and assigning a logarithmic non-uniformity to the step-sizes. This way, softer passages will be assigned more bits and finer resolution while louder passages will be assigned less (than uniform case) bits.
Futhermore, the distortion at the louder passages is perceptually less audible, due to psycho-acoustic principles, in particular to the masking effect; presence of a loud tone, suppresses the perception of nearby quiter tones, and makes them inaudible. Hence the quantization noise will be less audible while occuring with louder tones.


Estimating the Impulse Response of the Room Using Sweep Signal Microphone Recorded Signal (Input & Output of a Convolution)

this is the two-channel FFT method of spectrum analyzer:

$$ y[n] = h[n] \ \circledast \ x[n] $$

just make sure that the length of the FFT $N$ is at least as large as the length of sound $x[n]$ plus the expected length of the impulse response $h[n]$.  the length of sound $y[n]$ is also as long as the FFT.  you can round $N$ up to the nearest power of two.  just zero-pad everything to that length and then

$$H[k] = \frac{Y[k]}{X[k]} $$

is functionally true.

you might sometimes have to worry about division by zero, but if your driving signal $x[n]$ is sufficiently broad-banded (which a linear sweep or a maximum-length sequence is), then you don't have to worry too much about division by zero.

if you know $H[k]$, then you know $h[n]$.


1D filter for speed input and noisy position sensor

Sounds like an Extended Kalman Filter is the ideal candidate.

Bayes++ is a really simple drop-in library for implementing one, and the example here (https://github.com/Exadios/Bayes-/blob/master/PV/PV.cpp implements the question exactly (except it only takes in 1 observation)

Below is the code, in case it gets moved:

/*
 * Bayes++ the Bayesian Filtering Library
 * Copyright (c) 2002 Michael Stevens
 * See accompanying Bayes++.htm for terms and conditions of use.
 *
 * $Header: /cvsroot/bayesclasses/Bayes++/PV/PV.cpp,v 1.17.2.3 2005/04/07 06:39:22 mistevens Exp $
 */

/*
 * Example of using Bayesian Filter Class to solve a simple problem.
 *  The example implements a Position and Velocity Filter with a Position observation.
 *  The motion model is the so called IOU Integrated Ornstein-Uhlenbeck Process Ref[1]
 *    Velocity is Brownian with a trend towards zero proportional to the velocity
 *    Position is just Velocity integrated.
 *  This model has a well defined velocity and the mean squared speed is parameterised. Also
 *  the velocity correlation is parameterised.
 *
 * Two implementations are demonstrated
 *  1) A direct filter
 *  2) An indirect filter where the filter is preformed on error and state is estimated indirectly
 * Reference
 * [1] "Bayesian Multiple Target Tracking" Lawrence D Stone, Carl A Barlow, Thomas L Corwin
 */

#include "BayesFilter/UDFlt.hpp"
#include "BayesFilter/filters/indirect.hpp"
#include "Test/random.hpp"
#include <cmath>
#include <iostream>
#include <boost/numeric/ublas/io.hpp>

namespace
{
    using namespace Bayesian_filter;
    using namespace Bayesian_filter_matrix;

    // Choose Filtering Scheme to use
    typedef UD_scheme FilterScheme;

    // Square
    template <class scalar>
    inline scalar sqr(scalar x)
    {
        return x*x;
    }

    // Random numbers from Boost
    Bayesian_filter_test::Boost_random localRng;

    // Constant Dimensions
    const unsigned NX = 2;          // Filter State dimension   (Position, Velocity)

    // Filter Parameters
    // Prediction parameters for Integrated Ornstein-Uhlembeck Process
    const Float dt = 0.01;
    const Float V_NOISE = 0.1;  // Velocity noise, giving mean squared error bound
    const Float V_GAMMA = 1.;   // Velocity correlation, giving velocity change time constant
    // Filter's Initial state uncertainty: System state is unknown
    const Float i_P_NOISE = 1000.;
    const Float i_V_NOISE = 10.;
    // Noise on observing system state
    const Float OBS_INTERVAL = 0.10;
    const Float OBS_NOISE = 0.001;

}//namespace

/*
 * Prediction model
 * Linear state predict model
 */
class PVpredict : public Linear_predict_model
{
public:
    PVpredict();
};

PVpredict::PVpredict() : Linear_predict_model(NX, 1)
{
    // Position Velocity dependance
    const Float Fvv = exp(-dt*V_GAMMA);
    Fx(0,0) = 1.;
    Fx(0,1) = dt;
    Fx(1,0) = 0.;
    Fx(1,1) = Fvv;
    // Setup constant noise model: G is identity
    q[0] = dt*sqr((1-Fvv)*V_NOISE);
    G(0,0) = 0.;
    G(1,0) = 1.;
}

/*
 * Position Observation model
 * Linear observation is addative uncorrelated model
 */
class PVobserve : public Linrz_uncorrelated_observe_model
{
    mutable Vec z_pred;
public:
    PVobserve ();
    const Vec& h(const Vec& x) const
    {
        z_pred[0] = x[0];
        return z_pred;
    };
};

PVobserve::PVobserve () :
    Linrz_uncorrelated_observe_model(NX,1), z_pred(1)
{
    // Linear model
    Hx(0,0) = 1;
    Hx(0,1) = 0.;
    // Observation Noise variance
    Zv[0] = sqr(OBS_NOISE);
}

void initialise (Kalman_state_filter& kf, const Vec& initState)
/*
 * Initialise Kalman filter with an initial guess for the system state and fixed covariance
 */
{
    // Initialise state guess and covarince
    kf.X.clear();
    kf.X(0,0) = sqr(i_P_NOISE);
    kf.X(1,1) = sqr(i_V_NOISE);

    kf.init_kalman (initState, kf.X);
}

int main()
{
    // global setup
    std::cout.flags(std::ios::scientific); std::cout.precision(6);

    // Setup the test filters
    Vec x_true (NX);

    // True State to be observed
    x_true[0] = 1000.;  // Position
    x_true[1] = 1.0;    // Velocity

    std::cout << "Position Velocity" << std::endl;
    std::cout << "True Initial  " << x_true << std::endl;

    // Construct Prediction and Observation model and filter
    // Give the filter an initial guess of the system state
    PVpredict linearPredict;
    PVobserve linearObserve;
    Vec x_guess(NX);
    x_guess[0] = 900.;
    x_guess[1] = 1.5;
    std::cout << "Guess Initial " << x_guess << std::endl;

    // f1 Direct filter construct and initialize with initial state guess
    FilterScheme f1(NX,NX);
    initialise (f1, x_guess);

    // f2 Indirect filter construct and Initialize with initial state guess
    FilterScheme error_filter(NX,NX);
    Indirect_kalman_filter<FilterScheme> f2(error_filter);
    initialise (f2, x_guess);

    // Iterate the filter with test observations
    Vec u(1), z_true(1), z(1);
    Float time = 0.; Float obs_time = 0.;
    for (unsigned i = 0; i < 100; ++i)
    {
        // Predict true state using Normally distributed acceleration
        // This is a Guassian
        x_true = linearPredict.f(x_true);
        localRng.normal (u);        // normally distributed mean 0., stdDev for stationary IOU
        x_true[1] += u[0]* sqr(V_NOISE) / (2*V_GAMMA);

        // Predict filter with known pertubation
        f1.predict (linearPredict);
        f2.predict (linearPredict);
        time += dt;

        // Observation time
        if (obs_time <= time)
        {
            // True Observation
            z_true[0] = x_true[0];

            // Observation with addative noise
            localRng.normal (z, z_true[0], OBS_NOISE);  // normally distributed mean z_true[0], stdDev OBS_NOISE.

            // Filter observation
            f1.observe (linearObserve, z);
            f2.observe (linearObserve, z);

            obs_time += OBS_INTERVAL;
        }
    }

    // Update the filter to state and covariance are available
    f1.update ();
    f2.update ();

    // Print everything: filter state and covariance
    std::cout <<"True     " << x_true << std::endl;
    std::cout <<"Direct   " << f1.x << ',' << f1.X <<std::endl;
    std::cout <<"Indirect " << f2.x << ',' << f2.X << std::endl;;
    return 0;
}


Store a FFT with a minimal amount of data

In general, quantizing an FFT result is a lossy process.  Given that the FFT twiddle factors are transcendental functions, any finite storage size will result in adding some quantization noise; and that noise increases as the resulting format decreases in bit size.  For extremely sparse FFT results, it may be possible to compress those results, but those are statistically very rare cases.


Let me first state the commonly used definitions of minimum-phase and maximum-phase for discrete-time systems:

  A causal and stable discrete-time system is said to be strictly minimum-phase when all its zeros are inside the unit circle.

Of course, all its poles are inside the circle as well, otherwise it wouldn't be causal and stable.

  A causal and stable discrete-time system is said to be strictly maximum-phase when all its zeros are outside the unit circle.

Again, causality and stability requires that all the poles are inside the unit circle as well.

If you have zeros on the unit circle (and all others either inside or outside the circle), the system is usually referred to as loosely minimum-phase or loosely maximum-phase. A system with all its zeros on the unit circle is a degenerate case of that definition, in the sense that it is loosely minimum-phase and loosely maximum-phase at the same time. But in that case these definitions do not make much sense. These definitions are meant for systems that have some zeros inside (or outside) the circle and some on the circle.

Your question about the low pass frequency response can't be answered, because the concepts of minimum-phase and maximum-phase have absolutely nothing to do with the magnitude response of the filter. You can have a minimum-phase low pass filter, and a maximum-phase low pass filter with exactly the same magnitude response. The same is true for all other frequency selective filter characteristics (high pass, band pass, etc.).

Concerning your last question, it depends what you mean by shortest. If you mean the number of samples of the impulse response, then the answer is no. As an example, take a length $N$ FIR minimum-phase impulse response. If you flip the coefficients (i.e., time reverse them) you obtain a maximum-phase filter of the same length with the same magnitude response.

The minimum-phase property means that of all filters with the same magnitude response, the impulse response of the minimum-phase filter is maximally concentrated near time index $n=0$ in a mean square sense. I.e., if $h[n]$ is a minimum-phase impulse response, then

$$\sum_{n=0}^{N}|h[n]|^2\ge \sum_{n=0}^{N}|g[n]|^2\tag{1}$$

for all $N\ge 0$, and for any $g[n]$ with the same magnitude response as $h[n]$.


How to create matched "minimum phase" for a system of parallel FIRs?

The Hilbert "transform" or relation define one phase response for a given magnitude response, so you can't get both matched and minimum phase in your case.

Regards.


How to Calculate the Period of the Discrete Time Sequence from Taking Its DFT

Sample the signal into MATLAB.
Apply the fft function on it.

Plot the absolute value of the DFT signal (Its first half of samples).
Look for its peak, the index number of the peak multiplied by the sampling frequency divided by the number of samples will give the frequency.
Look at its inverse and your answer will pop.

If you upload a digital copy of the signal, I will show you in code.


You have 2 right-half plane zeroes : 0.012 and 18. The zero at +18 is "fast" and will not affect the performance much, but your slow zero at 0.012 will severely limit your performance. You can't cancel this right-half plane zero with a right-half plane pole in your controller, your controller output will be unbounded.
Is this homework or a real life problem?


Is it possible to use a Volterra series to generate subharmonics?

Check out this paper.  I would have made a comment but not high enough rep.

http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=1211087&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel5%2F81%2F27258%2F01211087

Looks like you need multiple in to get subharmonics in Volterra series

The abstract states "Subharmonic generation is a complex nonlinear phenomenon which can arise from nonlinear oscillations, bifurcation and chaos. It is well known that single-input-single-output Volterra series cannot currently be applied to model systems which exhibit subharmonics. A new modeling alternative is introduced in this paper which overcomes these restrictions by using local multiple input single output Volterra models. The generalized frequency-response functions can then be applied to interpret systems with subharmonics in the frequency domain."


More detail like PAK-9 asks, and example images with the pattern indicated, would help get better answers.  If you update your question with details and data, you're also more likely to get example implementations.

That said, this looks to be a candidate for a wavelet transform (like with pywavelets), which would allow for the stretching and scaling that you need to capture the variations of the pattern. Pre-computed wavelets can be browsed here, or you can roll your own (like the fast 2D CDF 9/7 lifting transform).

You suggested cross-correlation, which might be good for a first pass (see scikit-image template matching for ideas).  You might modify this to loop through windows to account for the stretching and scaling, but you'll want to normalize them if you intend to subsequently use a classifier.  Also, the relative width and spacing of the peaks seem conserved, which would lend itself to windowed supervised classification by many means, like SVM or Random Forests (scikit-learn has plenty of classifiers and examples for you).


Assuming $y$ and $u$ are scalars, then the first equation:

$$
u = -K \hat{x}
$$

is just a vector-vector multiplication and can be implemented in C as:

u = 0.0;
for (i=0; i<N; i++)
{
    u += xhat[i]*K[i];
}

assuming that $K$ is $1 \times N$ and that $\hat{x}$ is $N \times 1$.

The multiplication C*xhat can be implemented similarly (because I suspect that $C$ is really $1 \times N$ rather than $2 \times N$).

The multiplications L*(y-C*xhat) and B*u are vector-scalar multiplications which, for the latter, would look like:

Bu = 0.0;
for (i=0; i<N; i++)
{
    Bu += B[i]*u;
}

The multiplication  A*xhat  will be slightly different because $A$ is $N \times N$.


Differences between phase and group delay

Note that a constant group delay is not sufficient for a band-limited signal to exhibit no dispersion. It is the phase delay that needs to be constant. If the phase is affine, i.e., if we have

$$\phi(\omega)=a+b\omega,\qquad \omega>0\tag{1}$$

the group delay is constant

$$\tau_g(\omega)=-\frac{d\phi(\omega)}{d\omega}=-b\tag{2}$$

but the phase delay is not constant if $a\neq 0$:

$$\tau_p(\omega)=-\frac{\phi(\omega)}{\omega}=-\frac{a}{\omega}-b\tag{3}$$

And it is the phase delay that determines the delay of a sinusoid at the output of an LTI system. For $x(t)=\sin(\omega_0t)$, the output of an LTI system with magnitude response $M(\omega)$ and phase delay $\tau_p(\omega)$ is given by

$$y(t)=M(\omega_0)\sin(\omega_0(t-\tau_p(\omega_0)))\tag{4}$$

Only for a (strictly) linear phase (i.e., for $a=0$ in Eq. $(1)$) does the phase delay equal the group delay. In that case, the delay experienced by a sinusoid at the input of an LTI system becomes independent of frequency.


Convolution to do FIR with initial and final conditions!

Your benchmark is actually quite surprising! I would've expected filter to choose an appropriate implementation of convolution to filter the input data, but thinking about this, since filter needs to pass the state in a format that isn't dependent on the convolution implementation actually used, the whole process of generating the state might take longer than the actual filtering; note that this really is just speculation. (Assumption: filter with long input data uses FFT FIR implementation internally, but the state that filter returns are taps of a direct-form FIR implementation – you basically need to reverse quite a bit of the filtering you've done to get that)

  I already tried to overlap the chunks of data and mimic the "state" for the convolution but haven't been successful yet.

What you probably would like to have is a well-implemented FFT FIR! Since b doesn't seem to change, this would be a very efficient way of implementing the filtering operation.

The basic idea would be described as overlap-add filter (wikipedia), and the rough idea is this:

Pad your tap vector $b$ with zeros to length $N$. Get the $N$-point discrete Fourier transform of that using the FFT. Call the result $H$.
Take segments of length $L$ of the input signal, and pad them to $N$, transform as above.
point-wise multiply, since multiplication in discrete freq domain is equivalent to circular convolution in time domain.
inverse DFT (via IFFT)
By zero-padding enough, we now that a $L$ long part of the output signal can be generated by adding up the $N-L$ last samples of the last iteration to the start of the samples of this iteration.
save the $N-L$ last output samples for next iteration

personal addition: If you already know what filter you want, and you're actually after speed and don't care as much that you might need to leave Matlab to do the filtering (I personally cherish every moment that I don't have to use Matlab), look at this explanatory blog post that I've written, which actually does a lot of of band pass filters in parallel – you'd only need one of these, and you could choose a FFT FIR block to do the work; I've been doing 160 megasamples per second with these, on a "normal" PC. This is, as I'll frankly admit, based on GNU Radio and the VOLK acceleration library, because that's what I prefer for streaming samples DSP (which you seem to be doing in Matlab).


What is the bandwidth of a (real) sinusoidal tone, and pulse?

The spectrum of a continuous tone is, as you said, of the form $\delta(f-f_0) + \delta(f+f_0)$: 2 impulses at frequencies $f_0$ and $-f_0$.

As a lowpass signal, this is said to have bandwidth $f_0$ (the one-sided spectrum has components up to $f_0$).

As a bandpass signal, it has zero bandwidth (there's nothing around the carrier frequency $f_0$).

If you multiply the sine wave by a pulse, this makes it time-limited, and therefore frequency-unlimited. Infinite bandwidth in theory.

In practice, you must define some criteria for estimating your bandwidth. Examples are:

3 dB drop (of the sinc function around $f_0$)
10 dB drop
drop below the noise level


I need to see the code behind your graph because the numbers you put down are not making sense with that figure. Your graph should cover X values from 1-12501 if it's just FFT indices or 0-1250 for actual frequencies.

If you're looking to filter the two frequencies with the greatest magnitude no matter what the input signal is, you can find their frequencies by searching for local maxima. If you have the Signal Processing Toolbox, you can use this code to find the peaks and their indices, and then convert those indices into an actual radial frequency.

[pks,locs] = findpeaks(mag_fft_out);   % this fft should be 0-fs/2 ONLY
[throw_away,I1] = max(pks);            % find biggest peak and return index in peak vector
n1 = locs(I1);                         % find location of peak index in FFT vector
omega1 = 2.*pi.*n1./N;                 % N is the length of your original sequence in samples (T./fs)
pks(I1) = [];                          % remove primary peak from peak vector
locs(I1) = [];                         % remove primary peak from peak location vector
[throw_away,I2] = max(pks);            % find (remaining/second) biggest peak and return index
n2 = locs(I2);                         % find location of peak index in FFT vector
omega2 = 2.*pi.*n2./N;

This will also work for finding the peaks once or you could figure it out manually (Mk I Eyeball + Math)

After that you need to design a filter for each peak frequency. I highly recommend you use the MATLAB 1-D Digital Filter function filter(b,a,x) and choose your coefficients based on what characteristics you need. In all likelihood, a simple 2nd order bandpass will work for both. Look up how to design a digital filter. The article on Wikipedia is a great start and even has diagrams showing how the coefficients work.

I'm not sure this is what you're asking for. You say "delete" but it's not quite that simple. I supposed you could set all other values to zero except the peaks, and then if you want it back as a signal you inverse FFT it, but I've never actually tried that.


Smoothing signal / detecting bumps in a data stream

Effects of Averaging

Using a moving average filter will smooth out the irregularities in the signal. The noise becomes E/N where N is the length of the moving average filter. The side effect of using a MA is that the signal peaks become wider and shallower.

In addition, the frequency content of the signal will change. A moving average filter in the time domain is the same thing as convolving the frequency domain signal by a sinc function everything gets smudged out.

Peak Detection Algorithm
Peak detection is a common problem in 9/10 engineering problems. (not really, but a TON depend on them)

Typically this is what is done:

Median Thresholding

1) Look for all peaks in your signal. (i.e., a point that is larger than the two
   adjacent points
2) take this list of points and for each one of them compute:
   med_threshold = median(Peak,Width) + constantThresholmedian   where median is the
   median value of the data centered at "Peak" with Width being the number of
   points to look at.
       a) The Width(usually written as Lambda in literature) and constantThreshold
          (usually written as C) are determined by trial and error and using the ROC
          curve (Acronym below)
3) if the peak's magnitude is above this threshold accept it as a true peak.
   Else: Discard it, its a false peak
4) Generate a Receiver Operating Characteristic Curve(ROC) to how well the algorithm
   is performing.

Here is an example:

suppose we have the signal X = [ 0 0 0 0 1 3 **9** 2 1 1 **2** 1 1 ]
1) 9 and 2 are both potential peaks
2) Lets use a window of 5 and  a threshold =2
so at 9 we have [1 3 9 1 2] -> [1 1 2 3 9]  so Median(9,5) = 2
9 > 2 +2, therefor its a peak
Lets take a look at 2: [ 1 1 2 1 1] -> [1 1 1 1 2 ] Median(2,5) = 1
2 < 1+2, therefor it is NOT a peak.

Determining Frequency

Now that You've effectively found the time localization of the peak try to find their frequency:

1) Use the locations of the peaks to generate a pulse train
      a) this means create sum(Dirac_delta[t-L(n)]) where L(n) is the nth time that
      you've localized through median thresholding
2) Apply FFT Algorithm
3) Look for largest peak.

Alternate Frequency Estimation

1) Think of this like a beat in a piece of music (I learned about thresholding by
   researching Onset Detection.
2) Compute the average time distance between detected peaks.
3) now call your results BPM or PPM (pulses per minute)

Additional Avenues of Research

While you may be satisfied by the peaked signal as is, there are algorithms that are applied to a whole different beast of problems called Onset Detection.

Onset Detection is a big area in Music Information Retrieval Research. Its used to determine when a note is being played.

If you think of your tape head signal as a highly sampled signal, you can apply many of the algorithms that you would find in this paper:

http://www.elec.qmul.ac.uk/people/juan/Documents/Bello-TSAP-2005.pdf


what is the difference between the convolution and differentiation over image in image processing?

In this context, convolution refers to the operation of LTI filtering performed on the image by the filter impulse response. Filtering can have many purposes such as blurring, sharpening, noise reduction etc.

Certain applications requires that you compute an approximation to the derivative of the image data. This can be accomplished by filering the image with a specific filter kernel (impulse response) that's some sort of a high pass characteristics.

Therefore, computing the (approximate) derivative of an image can be accomplished by LTI filtering with a highpass impulse response, which refers to a convolution operation. This is the connection in between.


What does it mean that one frequency aliases another?

I assume that since you're watching lectures on sampling theory, then you have already grasped the core ideas of basic signal operations and Fourier transforms and signal spectrum plots.

Then, consider a signal $x(t)$ whose C.T. Fourier transform is $X(\omega)$. When this signal is multiplied by the (ideal) periodic impulse train $$ \sum_{k=-\infty}^{\infty} \delta(t-kT_s) $$ where $T_s$ is the sampling period (and therefore $F_s = 1/T_s$ is the sampling frequency in Hz.) the resulting signal is called as the sampled signal $$ s(t)=x(t) \sum_{k=-\infty}^{\infty} \delta(t-k T_s) $$

Then it can be shown (by applying properties of CTFT) that the spectrum $S(\omega)$ of the sampled signal $s(t)$ is the periodical extension of the spectrum $X(\omega)$ of the signal $x(t)$, i.e., $S(\omega)$ is obtained by shifting $X(\omega)$ to left and right in frequency by the amount $ k \times \omega_{s} = k \times 2 \pi F_s$.  ( Here I'm using radian frequency $\omega$ which is related to Hertz frequency $f$ by $\omega=2\pi f$ ) :
$$ S(\omega) = \frac{1}{T_s} \sum \limits_{k=-\infty}^{\infty} X(\omega -k 2\pi F_s) $$

For example consider the the shift for $k=-1$, then $X(\omega+2\pi F_s)$ is added on top of original spectrum $X(\omega)$. Now by setting $\omega=0$ in $S(\omega)$, it can be seen that $X(0+2\pi F_s)=X(2\pi Fs)$ has come to the origin in addition to $X(0)$. Setting $k=-2$ would show also that $X(4\pi F_s)$ has come to the origin. This way it can be deduced that those frequencies given by $\omega_k = k 2\pi F_s$ are all appearing at the origin in $S(\omega)$. This appearance of those frequencies $2 \pi k F_s$ at the origin (or more generally any frequency  $\omega_k = 2 \pi k F_s + w_0$ at the frequency $\omega_0$) is called as aliasing in the spectrum of $S(\omega)$


Is there a technical term for this simple method of smoothing out a signal?

What you've implemented is a single-pole lowpass filter, sometimes called a leaky integrator. Your signal has the difference equation:

$$
y[n] = 0.8 y[n-1] + 0.2 x[n]
$$

where $x[n]$ is the input (the unsmoothed bin value) and $y[n]$ is the smoothed bin value. This is a common way of implementing a simple, low-complexity lowpass filter. I've written about them several times before in previous answers; see [1] [2] [3].


What language or framework would you use to develop and test natively algorithms meant to be ported on an embedded hardware later?

I finally settled on python, and here is why:

I had no reason to limit myself to a specific language or framework. The target plateform was not specified, but I needed to have proof of concepts for the project, and quick.

I choose python because of its versatility, the profusion of libraries, the speed over other scientific calculation packages, and the scientific community.

And so far so good.


MATLAB Implementation of Karplus Strong algorithm with filter function?

The filter(b,a,x) function of MATLAB does not operate on symbolic variables. Your line of code

syms z

creates a symbolic variable z, in MATLAB environment, and later
you feed this variable into the filter() function through its second argument in this line:

H = filter(1,1-alpha*z^(-M),x);

where the filter(b,a,x) function, in this mode is expecting a numerical coefficient vector $a[k]$ associated with the DT system that's defined by the coefficient vectors $a[k]$ and $b[k]$.

In order to get the filter() function working, replace the symbolic expression

1-alpha*z^(-M)

with the numerical one, given integer M and real alpha

a = [1 zeros(1,M-1) -alpha]


Yes!

What is done in practice is simply disconnecting (or even better, terminating properly) the RF input of the receiver.

Then, you observe empiric noise variance separately on I and Q, and normalize so that they are equal.

This assumes that NF variations are benign, but it does work. If you can, you'd of course try to increase noise power.

USRP transceiver daughterboards can be RX IQ imbalance calibrated by letting TX transmit pseudowhite noise at a slight frequency offset (to average out potential TX imbalance) and observing the intra-board crosstalk. Same would work for any sufficiently random external signal. If you have a dominant external signal, intentionally tune a bit off, and average enough.


Discrete Wavelet Transform output: coefficients or FIR-filtered signals?

In one level DWT, each output of the low-pass or a high-pass can indeed be considered as signals. Thus each of those signals are subsampled by a factor of 2, and the same two-filter-subsampling is iterated on the low-pass output, several times (wavelet decomposition) at $L$ levels. Each final output of the different branches could still individually be considered as "signals", but they only make sense together, with respect to the input signal. Coefficients is a traditional name for them. Honestly, I don't know where it appeared originally.

Each represents information cast at a certain scale, and they don't really belong to the same domain. They are signals projected onto a nested wavelet subspace from a multiresolution analysis. They altogether belong to a global transformation. In other words: the values observed are signals with a heavy wavelet footprint (especially because of the subsampling). Analyzing them alone (without taking into account the wavelet projection) is theoretically dangerous, and even called "wavelet sin", or "wavelet crime", see for instance: Wavelet transform: How to compute the initial coefficients when only samples are available? (Strang and Nguyen, Wavelets and Filter banks, 1996, pages 232 sq.). However, people often may sin with practical efficiency, and for orthogonal non-redundant transformation, it often does not go too bad.

Anyway, due to the invertibility of the DWT scheme, you would have "more genuine" signals by: keeping some coefficients (setting the others to zero), especially from each "subband coefficient set alone", and performs the (unique) inverse DWT for each of the  $L+1$ sets of subband coefficients. From those $L+1$ inverted, you get $L+1$ reconstructed signals in the original domains, whose sum yields the original signal (linearity of the wavelet).

The "wavelet crime" often becomes more stringent when one uses union of wavelet bases, wavelet frames, etc. As for the MODWT, I would not claim (yet) for sure that I fully understand the concept. Its setting sis akin to undecimated, shift-invariant cycle-spinning or redundant discrete wavelets. While the DWT generates $N$ coefficients for a $N$-sample signal, whatever the levels, the MODWT generates about $(L+1)N$ coefficients. As far as I get it, MODWTMRA does a job similar to the one described above for the DWT: project the redundant wavelet coefficients back to th original domain, so that their sum yields the original signal back. Side note: for redundant transforms, the inverse is not unique.


transfer function of a sampler in the s domain

ADCs don't have a transfer function. An ideal ADC is an ideal sampler with a sampling period T.
Real ADCs are not ideal, they have some sampling delays, jitter, quantization noise, thermal noise, etc. Usually the sampling delays are really small compared to the dynamic of the systems, so we don't model them. Quantization noise and thermal noise should not be a factor unless the signal you measure is really small. And jitter is not an issue unless the signal frequency is really high.
So bottom line, if your ADC is really fast compared to the dynamics of your system and if the measured signal is inside the dynamic range of your ADC, you could pretend that the ADC is not there.


This is a really very nice applied signal processing question. Unfortunately your lack of DSP and physics expertise, builds some barriers on the way of your understanding. Hopefully this answer helps to tunnel through.

Now, your first question is about downsampling a discrete time signal. In fact there are already plenty of such questions here under the title of interpolators, downsampling etc. These are already answered with great detail. You better make a search on the site for more detailed answers but here once more I'll give a brief reply.

In downsampling a discrete-time signal (say from 360 hz to 60 hz in your case, by a facttor 6) the most critical thing to consider is the Spectral Utilisation of the original (360 Hz) signal. Based on Shannon-Nyquist sampling theorem criteria, this signal is assumed to contain meaningful (non-trivial) information in the band of frequencies 0 to 180 hz. However for many applications the frequency band of interest within the say 0-180 hz range will be a much smaller band, say for example from 0-20 hz. In such a case, the original signal is said to be "oversampled" (by a factor of 9 in this case). When a signal is oversampled by a factor of say M, then it is somewhat acceptable to discard every (M-1) samples and use the Mth one to represent the information content of the original signal without loss (except at the higher end of the band which gets to limit of aliasing). But this will/may produce some noise due to aliasing being on the very edge of it.

There won't be any aliasing, as long as the out-of-interest frequency band in the original signal contains no energy. In order to guarantee this condition of zero energy in the out of interest band, it is required to filter out those out of interest bands by a suitable filter of sufficient quality, before decimating the signal.

The quality of such a filter depends on, therefore, i-the spectral distribution of the original signal, ii-band of interest, iii-downsampling factor and other requirements such as real time processing or offline, accuracy, CPU power etc.

In your case, the answer lies in the following: the best downsampler should produce the most accurate representation of the band of interest signals and reject any irrelevant information for your steering hardware to process.

The short answer is this: as long as your system is configured to use 60 hz information (30 hz band) for the force-feedback processing, it seems that (and from my experience with mechanical acceleration signals, it is indeed so) this band of signals is adequate to represent average mechanical control. (as long as I can remember, most mechanical acceleration sensors do use at most a band of a few hundred hz). So just use a simple low-pass filter of cutoff frequency at pi/6 for better performance.

The second questions simply asks which is better to choose i-filtering 360hz into 20hz or ii-60 hz into 20 hz. The answer is the latter. It will be a much shorter filter to apply, and the transition band would also be easier.

When the latency of a few samples becomes so critical as in your example of real time forced-feedback control of a physical simulator, then high accuracy filters becomes impractical to use. May be you should completely discard using any filters at all, as you have pointed out the physical limitations of the wheel mechanics already prevents bulding of responses faster than 20 hz...


Does any paper use energy harvest to in the baseband?

No.

Your question neglects the fact that a mixer can't transfer the whole received power to baseband (mixers are typically active, and most of the baseband energy comes from the LO that the device has to generate itself).

You can use any nonlinear device, even one without an external power source, but the mixed output's energy will always be way lower than the energy of the inputs, simply by the pure math behind what mixing does.

In general, when asking such question, you'd start with "where does what energy come from?", and with a clean, and simplified, block diagram of the receive RF chain including the mixer, and that would have answered that. If you're researching power harvesting, this is the MAIN thing you'll need to understand!