I would follow a signal block diagram based solution for this problem.

First as suggested in the comments, it's very helpful to investigate a few values of $y[n]$ and $x[n]$ for some $n$ :

$$
\begin{align}
y[0] &= 0 ~~~,~~~ y[1] = x[1] \\
y[2] &= 0 ~~~,~~~ y[3] = x[2] \\
y[4] &= 0 ~~~,~~~ y[5] = x[3] \\
y[6] &= 0 ~~~,~~~ y[7] = x[4] \\
\end{align}
$$

with some experience, or by a direct investigation of this above sequence it can be seen that the system that produce $y[n]$ from $x[n]$ is the folowing :

$$ x[n] \longrightarrow \boxed{ \uparrow 2} \longrightarrow  w[n] \longrightarrow \boxed{ z^1} \longrightarrow y[n]$$

where the up arrow indicates an expansion by $2$ and the $z^1$ indicates a left shift (advance) by one sample.

Write down the DTFT relations between those signals $x[n]$,$w[n]$ and $y[n]$ :

$$
\begin{align}
W(\omega) &= X(2\omega) \\
Y(\omega) &= e^{j\omega} W(\omega) \\ \\
Y(\omega) &= e^{j\omega} X(2\omega) \\
\end{align}
$$

And relate $2N$ point DFT $Y[k]$ of $y[n]$ to $N$ point DFT $X[k]$ of $x[n]$.

$$
\begin{align}
Y[k] &= Y(\omega_k) = Y(\frac{2\pi}{2N} k) &, ~~ k=0,1,...,2N-1\\
Y[k] &= e^{j \frac{2\pi}{2N} k} X(\frac{2\pi}{N} k) &, ~~ k = 0,1,...,2N-1\\
Y[k] &= e^{j \frac{\pi}{N} k} X[k] &, ~~ k = 0,1,...,2N-1
\end{align}
$$

It can be seen that the length $2N$ DFT sequence $Y[k]$ consists of two fold replica of length $N$ DFT sequence $X[k]$ weighted by $e^{j \frac{\pi}{N} k}$. Note that $X[k]$ is periodic in $k$ by $N$, hence repeats twice during $k = 0,1,...,2N-1$.


Trying to identify the best method for IIR realizations

I don't think "regardless of application" works here. Every application has different requirements and hence the optimum choice would be different.

However in gross oversimplification

Direct Form I is the "safest choice". Direct Form II is slightly better since it has only half the amount of state variables. Stay away from Direct Form II and Transposed Form I since the transfer functions to the state variables are only determined by the poles and there can be massive gain scaling issues
Cascaded is the easiest choice. It's fairly simple to convert between poles and zeros and actual filter coefficients. Parallel requires calculating the roots of a large polynomial which is a numerically ill defined problem.
There are different ways to pair poles and zeros to create second order sections. Best one is this : http://scipy.github.io/devdocs/generated/scipy.signal.zpk2sos.html
Use floating point whenever possible. Fixed point IIR requires very careful analysis and application specific implementation.
Whether optimization is worth it really depends on the application. Requirements can be all over the place: precision, speed, memory, mips, ease of us, stability against outlier condition, etc.


Relating the DFT to the FT (CTFT) is a big issue.  Let's start with the basic definitions without any domain specification using a $\frac{1}{N}$ normalization on the DFT.
$$ FT(x(t))(f) = \int x(t) e^{-i2\pi t f } dt  $$
$$ DFT(x[n])(k) = \frac{1}{N} \sum x[n] e^{-i 2\pi \frac{n}{N} k } $$
They are obviously similar.  Assuming they cover the same interval, we can map the domain variables to each other.
$$
\begin{aligned}
t &= \frac{n}{N} = \frac{1}{N} \cdot n \\
f &= k  \\
( \Delta t  &=\frac{1}{N} ) \to dt
\end{aligned}
$$
[Edit: Does anybody want to argue the (1/N) doesn't belong in the DFT definition?]
The units of $f$ and $k$ are cycles per the interval.
Of course, we are assuming the sequence of values matches the function at every point.
$$ x[n] = x(t) $$
The FT is defined on a domian of $-\infty$ to $\infty$, and the the DFT is defined on a domain of $0$ to $N-1$.  The corresponding $t$ values are $0$ and $\frac{N-1}{N}$.  We can now include the domains in the definitions for a direct comparison on the same interval.
$$ FT_{DFT}(x(t))(f) = \int_{0}^{\frac{N-1}{N}} x(t) e^{-i2\pi t f } dt  $$
$$ DFT(x[n])(k) = \frac{1}{N} \sum_{n=0}^{N-1} x[n] e^{-i 2\pi \frac{n}{N} k } $$
To make the $FT_{DFT}$ transform fit the $FT$ transform, we have to introduce a rectangular window function on the interval.
$$ w(t) =
\begin{cases}
= 0 & t < 0 \\
= 1 & 0 \le t < 1 \\
= 0 & 1 \le t \\
\end{cases}
$$
$$ FT_{DFT}(x(t))(f) = \int_{-\infty}^{\infty} w(t) x(t) e^{-i2\pi t f } dt  $$
Please notice that there is no window function in the DFT! The window function is in the FT.
From these definitions, it is very clear that the FT is the limit of the DFT with a $\frac{1}{N}$ normalization factor as $N$ goes to infinity.
$$
\begin{aligned}
\lim_{N \to \infty} DFT(x[n]) &= FT_{DFT}(x(t)) \\
\lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} x[n] e^{-i 2\pi \frac{n}{N} k } &= \int_{-\infty}^{\infty} w(t) x(t) e^{-i2\pi t f } dt \\
\end{aligned}
$$
This is the "right" way to view it.  An integral is the limit of a summation.
It is possible to to view the DFT as a sampled version of the FT by using a second window "function" (loosely speaking) of a train of Dirac deltas.  But this is not the definition of the DFT, requires advanced Real Analysis understanding, and you are merely getting back to the definition by working backwards.
Some EE professors neglect to mention this, some just don't know it.
So, just like you can approximate an intergral with a summation, you can approximate an FT on a finite interval with a $\frac{1}{N}$ normalized DFT.
To get intervals other than $[0,1)$ on $t$, call your new time $\tau$ and rescale it.
$$ t = \frac{\tau}{T} = \frac{n}{N} $$
Feel free to shift it too.
$$ t = \frac{\tau- \tau_0}{T} = \frac{n}{N} $$
This makes the sequence definition equivalent to
$$ x[n] = x_{\tau}(\tau) = x_{\tau}\left( \tau_0 + \frac{T}{N}n \right) $$
The DFT does not care what your sampling rate is. aka $\frac{N}{T}$ with these names nor where you place the frame, $\tau_0$.

Followup:
Questions from the OP's comment under the post question:
(All answers assume you are dealing with a pure real tone with a slighlty varying envelope.)
"How do you get rid of the influence of the number of samples and the range over which y(x) is defined? "
By framing your interval on a whole number of cycles.  This is the principal justification in my mind for a $\frac{1}{N}$ normalization factor of the DFT.  On a whole number of cycles (for a pure steady tone) a single conjugate pair of bins are non-zero.  The value in the bin is independent of N when normalized.  Thus the reading you get matches the continuous case exactly and is independent of the sample count.  Also, there is no need to calculate the whole DFT (take a FFT), a single bin calculation will do.  With a consistent frequency having a set of basis vectors is most efficient.  With a varying frequency, a Goertzel calculation yields the equivalent more efficiently.
"How do you obtain the correct amplitude of the Fourier transform of the function y(x)? "
With a single bin, the magnitude of the bin is half the amplitude of the tone.  A small DFT frame size of 2 to 3 cycles in length ensures that the reading is fairly accurate on that interval.  The reading will be that of the best fit sinusoid on that frame.
"How do you properly take into account the number of samples and the range over which y(x) is defined? "
Not sure what you mean by this.  The sample rate and signal frequency define how many samples per cycle you will have.  Decide the number of cycles you want, multiply that by the number of samples in a cycle and you get the sample count (N) for that frame.  If you are calculating a single bin, there is no special advantage to powers of two or any other factor.


Find a sinusoidal fit from data with fft

Depending on kind of signal you are dealing with I can think of the easiest method to do so.

Take the DFT of your signal (remember about proper scaling of values).
Find the dominant frequency $f_0$ and its amplitude $A$ at magnitude spectrum.
From phase spectrum find the corresponding phase shift $\theta$ for this given frequency - usually it is given in radians.
Reconstruct your signal as: $x(t)=A\sin(2\pi f_0t+\theta) $

Some quick and dirty code is included below. No warranty of working whatsoever.

#!/usr/bin/env python

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

if __name__ == "__main__":
    # Signal and analysis parameters
    fs = 4000.0    # Sampling frequency
    T  = 2         # Duration of a signal
    f0 = 100       # Fundamental frequency
    A  = 1         # Sine amplitude
    #th = np.pi/4   # Phase shift in radians (45 degrees)
    th = 0

    # Time vector
    t = np.arange(0, T*fs)/fs

    # Generate the signal with higher harmonics
    x = A*np.sin(2*np.pi*f0*t + th) + \
        A*np.sin(2*np.pi*2*f0*t + 2*th)/8 + \
        A*np.sin(2*np.pi*3*f0*t + 3*th)/16

    # Add some noise
    x += 0.5*np.random.randn(x.size)

    # Number of all samples
    N = x.size

    # Apply the window
    win = sig.gaussian(N, N)
    x = x*win

    # Perform the DFT
    X = np.fft.fft(x)
    # Get the frequency vector
    freq = np.linspace(0, fs-1/fs, N)

    # Get the phase and magnitude
    Xmag = np.abs(X)/np.sum(win)*2
    Xph  = np.angle(X)

    # Magnitude in logarithmic scale
    Xmagl = 10*np.log10(Xmag/Xmag.max())

    # Find the amplitude of the dominant frequency
    ind_max = np.argmax(Xmag[0:np.ceil(N/2)])

    # Extract the amplitude ...
    A_ = Xmag[ind_max]
    # ... frequency ..
    f0_ = freq[ind_max]
    # ... phase
    th_ = Xph[ind_max] + np.pi/2

    # Reconstruct the signal
    x_ = A_ * np.sin(2*np.pi*f0_*t + th_)

    # Do some plotting
    plt.subplot(311)
    plt.plot(t, x)
    plt.grid(True)
    plt.hold(True)
    plt.plot(t, x_, 'r')
    plt.title('Time domain signal')

    plt.subplot(312)
    plt.plot(freq[0:N/2], Xmagl[0:N/2])
    plt.grid(True)
    plt.title('Magnitude')

    plt.subplot(313)
    plt.plot(freq[0:N/2], Xph[0:N/2])
    plt.grid(True)
    plt.title('Phase')
    plt.show()


How to determine poles and zeros of the z-transform?

Just multiply numerator and denominator by $z^2$ to obtain

$$H(z)=\frac{z^2}{(z-\frac12)(z-2)}$$

from which you see that there's a double zero at $z=0$. Note that the number of zeros and poles is always equal if you include poles and zeros at infinity.


Connecting edges detected by an edge detector

This might not be complete solution, but will give you good direction.

Basically, what is the key criteria of to say that edges match? That "locally" the gradient of the edge matches and to some extent the distances are reasonable against how long the edge is continuous.

If you have geometric edges, like long straight lines, Hough will do very seamless job right away. But this doesn't work when edges are arbitrary curve. In this case, you can still think of curve as roughly piece wise constant segments (good enough in your case), hence, you should take Hough locally. i.e. you can take a small portion of image (say a block) compute the Hough and identify some peaks. Based on this, you can identify that de-houghing is creating meaningful gap, if so, keep it or move on.

Once smaller gaps are filled up, you can extend the same to take up on a larger size. the peaks will be more, but you may select fewer from it.


Why do Octave and MATLAB both automatically shift the gray values in this image?

Your image is an indexed image. Meaning it takes integer values that are supposed to be mapped to colours on screen via a particular colourmap.

You can see this if you read your image in as:

[I, M] = imread('dots_3gray.png');

where I becomes your 'indexed' image (which seems to contain indices from 0 to 158), and M is the colormap used to interpret it.

If you look at the documentation of imshow you will see that for indexed images specifically, you can specify which colourmap to use; in this case M to view it as expected:

imshow(I, M)

If you don't want this to be an indexed image, and you'd prefer to convert it to a standard grayscale image (i.e. in the range [0,1]) so that you can process it further, you can perform this conversion using the ind2gray command:

G = ind2gray (I,M); % convert to uint8 grayscale image in the range [0,255]
G = im2double (G);  % convert to double precision image in the range [0,1]

See my answer to a similar question in stackoverflow for more info.


Frequency is more likely decreasing (period is getting longer). This is a sweep sine signal or as some people used to call it - chirp tone. No point in rewriting the literature, so here are some links with very good explanation for both linear and exponential frequency change.

For linear sweep you get:

$$s(t)=\sin\left(\omega_1t+\frac{\omega_2-\omega_1}{T}\cdot \frac{t^2}{2}\right) $$

where:

$\omega_1 $, $\omega_2$ - initial and final angular frequency respectively

$T$  - sweep duration.

  Meng Q., et al. - Impulse Response Measurement With Sine Sweeps and
  Amplitude Modulation
  Schemes

  Farina A. - Advancements in impulse response measurements by sine
  sweeps


How to derive the answer to this convolution problem?

HINT: what is the highest power of $x$ after multiplying $A(x)$ and $B(x)$? This gives you directly the value of $N_3$. Then rewrite $B(x)$ as

$$B(x)=\sum_{n=0}^{2N_2}\hat{b}_nx^n$$

and you can use normal convolution of $a(n)$ and $\hat{b}(n)$ to derive $c(n)$. Now you just need to express $\hat{b}(n)$ in terms of $b(n)$.

EDIT: OK, so here's the solution:

We have

$$C(x)=\sum_{n=0}^{N_3}c_nx^n$$

Since $C(x)=A(x)B(x)$, the highest power of $C(x)$ must be $N_3=N_1+2N_2$. You can rewrite $B(x)$ as

$$B(x)=\sum_{n=0}^{2N_2}\hat{b}_nx^n$$

with

$$\hat{b}_n=\begin{cases}b_{n/2},&n \text{ even}\\
0,&n \text{ odd}\end{cases}$$

Now the coefficients $c_n$ can be written as the convolution of $a_n$ and  $\hat{b}_n$:

$$c_n=\sum_k\hat{b}_ka_{n-k}$$

For the index $k$ in the above sum we have the following constraints:

$$0\le k\le 2N_2\quad\text{and}\quad 0\le n-k\le N_1$$

which results in the summation limits

$$c_n=\sum_{k=\max\{0,n-N_1\}}^{\min\{n,2N_2\}}\hat{b}_ka_{n-k},\quad 0\le n\le N_3$$

If you sum only over even $k$, you can replace $\hat{b}_k$ by the coefficients $b_{k/2}$:

$$c_n=\sum_{k=\max\{0,n-N_1\},k\text{ even}}^{\min\{n,2N_2\}}b_{k/2}a_{n-k},\quad 0\le n\le N_3$$

which can be rewritten once more as

$$c_n=\sum_{k=\lceil{\max\{0,n-N_1\}/2}\rceil}^{\lfloor{\min\{n,2N_2\}/2}\rfloor}b_{k}a_{n-2k},\quad 0\le n\le N_3$$


how to do this low pass filter by hand on paper?

You can use very simple low pass filter - Moving average. Filter is quite simple - it is sum of N sample of signal divided by N. If length of filter will be equal to   period of carrier signal then performance of this filter will be good.


What's the relation between frequency band of $X(j\omega)$ and $\Phi_{xx}(j\omega)$?

The autocorrelation of $x(t)$ is
$$r_x(t)=x(t)\star x(-t)\tag{1}$$
where $\star$ denotes convolution. Taking the Fourier transform of $(1)$ gives
$$S_x(\omega)=X(\omega)X^*(\omega)=|X(\omega)|^2\tag{2}$$
$S_x(\omega)$ is the energy density of $x(t)$, and according to $(2)$ it equals the squared magnitude of the Fourier transform of $x(t)$. So if $x(t)$ is band-limited, both $X(\omega)$ and $S_x(\omega)$ are zero outside the signal's bandwidth.
Note that a deterministic continuous signal which has a Fourier transform (represented by an ordinary function) is usually an energy signal, which doesn't have a power spectrum (only an energy density).


The paper referenced in your link seems to be this one.

Of particular interest there is Table 1 (included below). The accuracy rates aren't great, though they are better than other approaches.


How to measure the "correctness" of a pink noise generator?

Simply observe the noise for a longish while and estimate a PSD of it – for example, simply by doing an FFT and observing the magnitude of that, and calculating the mean square error to the theoretical (triangular) PSD of pink noise.

That can be easily implemented only using Python/numpy (fft, abs, mean are all implemented in numpy).

Another, pretty sensible, method would be designing a filter bank that exploits the fact that octaves always contain the same amount of power in pink noise (which is why you'd use it in e.g. audio room sounding).

scipy.signal contains ample FIR filter design methods suitable for designing filters.

The band from 10-100 Hz should contain the same power (i.e. average sample square) as 100-1000 Hz, and 1000-10000 Hz, for example. Make sure scale the transition widths along with the passband widths. You don't have to go factors of 10, you can e.g. use 10-20 Hz, 20-40 Hz, 40-80 Hz … and verify that they all these contain the same power.


What's the best technique to determine if the contents of two images are the same?

A simpler approach (at least to start with) would be to evaluate the absolute difference between the two images and specifically, its mean value or maybe even just the sum of it (across all pixels), and then take its histogram across your sample of pairs of images.

The histogram will follow a bimodal (or multimodal) distribution with two main modes: The lower mode corresponding to the pairs of images with slight differences and the higher mode (higher values of difference, therefore more error) to pairs of images with drastic differences. You can then use this information obtained from the histogram to work out a threshold. For more information please see: http://en.wikipedia.org/wiki/Multimodal_distribution

A few details:

You don't have to work on the colour images for this, you can simply
convert them to greyscale. The way you do this will depend on the
tools you are using.
The difference image that you depict in your post seems to be a
bitmap, the absolute difference between two grayscale images would
also "look like" a grayscale image, please see:
http://homepages.inf.ed.ac.uk/rbf/HIPR2/pixsub.htm

Hope this helps, happy to amend the response if its not exactly you are looking for, I was a bit surprised to see you going for edge detection straight from the beginning so maybe I have not understood exactly the definition of "frame difference" you might be thinking about.


Can anyone explain me this MDFT polyphase filter bank?

It's input agnostic, anything will work just as it would with any other real valued prototype filter. I've implemented polyphase filter banks this on radar systems in practice, where we're operating on complex data, both pulse compressed and uncompressed. Filter banks like these have loads of applications due to the i inherent design and theoretical speed.
Polyphase filter banks typically create a phase shifted "copy" of the original version to achieve better reconstruction. From just scanning the code, it looks like this is what he's doing with his X1 and X2 vectors. You'll notice if you step through the code and check out the spectrum, the sub-bands will have no passband overlap. It's even more obvious if you look at the shifted versions of the filter frequency response. Since the "copy" of the original signal was phase shifted, it will also have non-overlapping segments within itself, but tougher the original and copy combine to fully cover from -pi to pic on a normalized frequency scale.

Typically the processing chain is this: get input, create a copy, phase shift the copy (everything for here on out is done twice, one for the original and one for the phase shifted copy; let's call these upper (original) and lower (copy)), apply prototype polyphase filter to the upper and lower signals, apply a DFT to do the polyphase magic, and now you're all channelized. Synthesis is pretty much just this in reverse.

The part people struggle with is typically the polyphase filters and the use of the upper and lower signals, which seems to be the case for this question . The filter bank itself is not complicated but the math behind the polyphase can be if you don't have a multi-rate signal processing background . The prototype filter is typically designed to make say 32 sub-bands, but since we have are upper and lower signals which each have 32 sub-bands, we really have 64. I'll leave the math behind why polyphase filters work to a textbook since that wasn't in your question.

See above, has to do with the phase shift.  Also go a few lines up and you'll see a negation. Check out those indices. The subtraction you're seeing may not be doing exactly what you think it's doing. Try stepping through the code and checking it out.


Usually, people tend to go with the second option, since its a joint approach. In either case, the methods fusing these two data modalities (RGB+Depth etc) together are termed multi-modal methods.

Also, another 3rd option is not to fuse on the feature level but to jointly minimize an energy function or jointly voting for the correct hypothesis.

One successful algorithm which works on RGB+Depth is LineMOD. The publications are here and here. (Modified Line). LineMOD takes the second path, and fuses the depth and RGB features. It also has an OpenCV implementation (a basic one).


To perform carrier recovery for FSK, if the data is random and can assumed to be equiprobable, then you can take the mean of the derivative of the phase versus time to determine the carrier offset.

Optionally the phase change from one sample to the next can be estimated using using a complex conjugate product of successive samples (for small angles, the imaginary result of the complex conjugate product of two normalized vectors is approximately the phase). Since this is the change in phase over a time interval, this result will be directly proportional to frequency:

$$f \propto Im(s(n)s(n-1)^*)$$

Where s(n) is the normalized complex FSK time domain signal.

Note in terms of I and Q samples this is simply:

$$I_{n-1}Q_n-I_nQ_{n-1}$$

The average of this result is proportional to the frequency of
the carrier and the deviation from average is the demodulated FSK signal. (The carrier component assuming a fixed offset can therefore be removed by using a simple DC Nulling Filter such as the one shown in this post: What does correcting IQ do?)

Since a small angle approximation was used, this provides a linear demodulation when the sampling rate is sufficiently high such that the phase change between two samples is a small angle such that $sin(\theta) \approx \theta$. For FSK that linearity is generally not as critical so would also work with much larger phase angles.


Best window before FFT ? (for a signal consisting of 2 tones used for phase measurement)

Using an FFT to measure phase for just two tones results in a lot more processing that doing the following alternate approach that can be either streamed or processed in blocks. No windowing is needed:
Apply the received signal as the input to two multipliers. Apply a normalized local copy of one tone as the second input to one of the multipliers and a normalized local copy of the other tone to the other. The low pass filter of each output will be an estimate of the phase given by$A\cos(\theta)$.
For full 360 degree phase resolution, use a complex local tone and two complex multipliers with two outputs (I, Q) with the phase determined using atan2(Q, I).   (This is essentially the operation in two bins of the DFT).


Simulate a simple multipath channel for OFDM Modulation test. MATLAB

Wav files expect values from -1 to +1.  You are probably exceeding those limits when you filter the signal through the channel model.  After the channel filter add the following:

after_channel = after_channel / max(abs(after_channel));

I suspect, though, that you will still have a problem.  I believe that wavwrite ignores the imaginary part of complex data, so you will be losing that part of the signal.


How to normalize output of $\tt scipy.signal.correlate$

When you say normalized cross-correlation I guess you mean the Pearson correlation.
Anyways you just divide the cross correlation by the multiplication of the std(standard deviation) of both signal, or more conveniently:
$ \rho_{xy}  =\frac{<x,y>}{\sigma_x\sigma_y}$

and in code:

x1 = x1/x1.std()
x2 = x2/x2.std() and then as you did it


DFT N-point circular convolution of $\cos(2 \pi n/N)$ and $\sin(2 \pi n/N)$

HINT : $\cos(2\pi n/N) = 0.5e^{j2\pi n/N}+0.5e^{-j2\pi n/N}$. This itself is a IDFT equation, where $X_1[k] = 0.5\delta[k-1] + 0.5\delta[k+1]$.

Similarly for $\sin(2\pi n/N)$, $X_2[k] = (0.5/j)\delta[k-1]-(0.5/j)\delta[k+1]$.

So $Y[k] = (0.25/j)\delta[k-1] + (0.25/j)\delta[k+1]$, and $y[n] = (1/2) \times (0.5/j)e^{j2\pi n/N} -(1/2) \times (0.5/j)e^{-j2\pi n/N} = (1/2) \sin(2\pi n/N)$

To expand mathematically
$$
X_1[k] = \sum (0.5e^{j2\pi n/N}+0.5e^{-j2\pi n/N})e^{-j2\pi kn/N}$$ $$= \sum (0.5e^{j2\pi (1-k)n/N }) +(0.5e^{j2\pi (-1-k)n/N })$$ $$= 0.5\frac{(e^{j2\pi (1-k)N/N} - 1)}{(e^{j2\pi (1-k)/N} - 1)} + 0.5\frac{(e^{j2\pi (-1-k)N/N} - 1)}{(e^{j2\pi (-1-k)/N} - 1)}$$ $$= 0.5\frac{e^{j\pi(1-k)}\sin(\pi(1-k))}{(e^{j\pi/N(1-k)}\sin(\pi/N(1-k))} + 0.5\frac{e^{j\pi(-1-k)}\sin(\pi(-1-k))}{(e^{j\pi/N(-1-k)}\sin(\pi/N(-1-k))}
$$
You can use similar steps to compute $X_2[k]$. You need NOT proceed to multiply from here. Look at the numerator term of each fraction. $\sin(\pi(m-k))/\sin(\pi/N(m-k))$ is zero when $k \ne m$. For $k=m$, even before you reach the geometric series sum, you can see all terms in summation are 1. So $X[k] = N$ when $k=m$. So $X_1[k] = N\times 0.5\delta[k-1] + N \times 0.5\delta[k+1]$ . Though I am bit confused about the $N$ factor.


Separate Simple Signal into Exponential Components

Your problem seems to be finding the analytic signal to a given real signal.

What you apparently need is the Hilbert transform. This will allow you to "recreate" the imaginary part of $y(t)$ from the real part alone, which is essentially your measurement $x(t)$.

Note that the transform can only be approximated as it is a acausal filter, and only works best when the signal $x(t)$ is available in the computer for all $t$.

A simple implementation is to Fourier transform $x(t)$ and set the negative frequency components to zero. After transforming back you obtain (the complex signal) $y(t)$.


Multiple short-range impulse responses aggregated to simulate a long-distance impulse response

if I were to convolve a source signal twice, would that result in a reasonable approximation of what the sound would be at 100m distance?

Yes. For small number of cascades this is a reasonable assumption.
However, you can't cascade the 10 meter response 10 times to get a single 100 meter response. While the general shape and structure of the response from 20m to 30m is the same as that from 30m to 40m, the phase spectra and the fine structures of the reflections are uncorrelatated. If you cascade the same impulse response over and over again, you get spectral thinning and it's starting sound "metallic" and "hollow".
You can try to decorrelate short repsonses from the original or just to start with a parametric model that has tree density, tree reflectivity, leaf/brush absorption and obvioulsy distance as tweakable parameters.


Whether Poisson noise will be dominant source of noise in low light conditions?

The Hasinoff paper says:

  For small photon counts, photon noise is generally dominated by other
  signal-independent sources of noise.

From my own experience with cameras I know that that is true. Even if the lens cap is on, the pixels will have dark current noise that is both spatially patterned and to a lesser extent varies by time. See for example this Nikon D300 dark frame:

The Alter et al. paper only considers the low-light case. They do not take into account any time-varying component of the dark current noise. They do mention an "offset voltage", and take it into account in their noise model. Maybe their sensors are just very good with an almost perfectly constant offset voltage for each sensor element.

The dark current noise is not photon noise, because by definition there are no photons, but it is still Poisson-distributed (see Justin Charles Dunlap, Characterization and Modeling of Nonlinear Dark Current in Digital Imagers, 2014), independent of the signal. Because there are these two types of noise sources with Poisson statistics, I don't think it is wise to call photon noise the Poisson noise like Hasinoff does.


In a typical quantization like that seen in an ADC the output would saturate any value > 3 to 3 (also known as clipping). Note that if you're considering a software quantization, it's also possible that some quantizers (although IMO poorly designed ones) could cause out of range values to wrap around and switch signs.


Why every musical instrument sounds different than any other while listening?

Actually there are classes of musical instruments whose members will have similar sounds such as the string-class of violin, viola, cello or the class of woodwind instruments with flute,picolo,pan flute etc. or sets of drums or the class of percussions including a grand piano, stage piano, console piano etc...

So its clear that each class has a very distinctive sound as an average person can easily classify them while within a given class it is more difficult for an average person to identify every member (but a person with a musical ear will be able to distinguish almost every type of instrument from each other, thats one of the primary differences between a musically able person from a not so one)

Coming to what makes all those instruments (both the classes and members)  sound different from each other is their physical structure and vibration mechanisms. These are specific to almost every class and to each member based on their material, geometry and driving mechanisms so that the resulting acoustical waves will have specific harmonic characteristics for each instrument (indeed this is true for any object of the nature) which will lead to distinguishable timbre / tone differences by humans.

Coming to the ability of simulating the acoustical charactheristics of an instrument by another one. There are already electronic analog/digital synthesizers which can mimic many of the instruments to some degree. There are also more sophisticated acoustical modeling programs which offer the complete re-creation of the acoustical properties of any imaginable or existing object/instrument. I'm not experienced in their results so you can look further.


Having images of a signal is very different from having a frequency offset in the signal. Having an image means the signal spectrum is replicated at two or more places.Having a frequency offset in a signal means the spectrum is just shifted in frequency.

First let us look at it from the point of view of the innocent receiver who does not know that such a catastrophe has occurred that the offset itself is equal to the channel frequency. Simply speaking the systems would break down with IQ phase reversed.

If however we have realistic frequency offsets like that of Doppler shifts which for typical cellular frequencies in current wireless is of the order of few hundred Hz then this is usually measured at the receiver and the clock is adjusted. This offset is usually measured by correlating time domain samples after ADC in baseband, since a linear shift in frequency domain equals phase shift in time domain.


this paper was done long before i had MATLAB.  the drawings are poor, but the math (at least in this revision, which is what you should use) is spot on.

send me an email address (to my audioimagination.com, not the wavemechanics.com on the paper) and i will send you a very short C file that shows how to generate and crossfade the wavetables in the synthesis process.

about analyzing a sample of a note (like from a wav file), you need to have a good sliding pitch-detector and a good interpolator (like maybe a good Kaiser-windowed $\operatorname{sinc}(\cdot)$-based interpolator) to extract single-cycle waveforms out of the note and put them into wavetables.

or you can create wavetables in memory from the POV of additive synthesis, simply by adding up the sinusoids, but it might be easier with the inverse FFT (and this is when knowledge of the periodicity of the DFT can be used to advantage).  i think your life will be easiest that the active wavetables in memory that are being used in the real-time synthesizer, those wavetables should all be the same power of 2 (and a nice big power of 2, like 1024 or 2048 or 4096).  notes that are higher in pitch will have fewer harmonics than notes lower in pitch, but will have the same number of samples per wavetable.  (and if the number of samples greatly exceeds twice the number of harmonics, then you can get away with linear interpolation during synthesis.)

wavetables that are being stored and not active, need not have so many samples.  as long as the number of samples in the periodic wavetable is more than twice the highest harmonic number, then, at least theoretically, it is enough.  but when a wavetable is made "active" and loaded into memory, it should be expanded and interpolated out to be whatever the large power of 2 is.

oh, and Nigel Redmon has code somewhere.


Architecture for high sample-rate FIR/IIR filter on FPGA or ASIC?

Short answer: If your sampling rate is higher than your throughput, you must use some kind of polyphase decomposition of your filters to work. So if you can't do that, you can't use a 100 MHz-clockable design to filter something with a sample rate of 1000 MS/s. It's as simple as that.

However, your example is exactly the standard example of a filter that can be implemented as a polyphase decimator: you get in 1000 MS/s, you filter to 100 MHz of that, and decimate by ten (note: assuming complex sampling here). You can do that, because your bandpass restricts the bandwidth that passes through to $\frac1{10}$ of the input bandwidth.

Anyway, even so, this will very probably not work out – there's simply zero margin for "real world" effects in your 100MHz clock to produce 100 MS/s, so you should really look into a design that can work faster.

Also, it's usually not that you buy an FPGA and can say "this works at 100MHz, no matter what I do on it", but you know its maximum clocking (which typically is much higher than 100MHz for devices that can take in 1GS/s via SERDES these days), and the maximum clock rate achievable depends on the DSP you do – so you can typically just ask your FPGA design tool to generate a FIR that is optimized for high rate (rather than e.g. for low power or small ressource consumption), and you'll get a different possible rate depending on what FIR architecture you chose.


Maximum likelihood estimation complexity computation

Ordinary Least Squares problem
Your
$$\hat x = \arg\min_x \sum_{n=1}^{N_r} \left\lvert y_n - h_n x\right\rvert^2$$
is just a way of saying
$$\hat x = \arg\min_x \left\| y - Hx \right \|$$
and that is called an Ordinary Least Squares estimation problem.
Its objective function is of a quadratic form, so we don't have to actually iteratively optimize or anything. There's an analytic solution to the above equation:
$$ \hat x = \left(H^T H\right)^{-1}H^T y\text.$$
As you can see, this has absolutely nothing to do with the constellation order so far – and it doesn't have to. The constellation order only makes a difference to the decision complexity, and by the orthogonality principle, that's a step to be done after you've estimated the transmit signal.
Don't invert that matrix
Now, it's practically never desirable to calculate a matrix inverse, so the thing that you actually solve is
$$\left(H^T H\right)\hat x = H^T y\text.$$
How complex the solution to that is depends on the properties of your MIMO channel matrix!
The quadratic / full-rank MIMO case
In the (very special) case that $H$ is normal (that implies among other things quadratic, ie. as many receive as transmit antennas $N_r=N_t=:N$), you do that by applying a QR decomposition to $H^T H$:
\begin{align}
\left(H^T H\right)\hat x &= H^T y\\
QR \hat x &= H^T y\\
R\hat x &= Q^T H^T y\tag{1}\\
&= \nu y\tag{2}
\end{align}
Since $R$ is a triangular matrix, $(1)$ s the point where you just backsubstitute. Luckily, $R$, $Q^T$ and $H^T$ only need to be calculated once ever, so we can estimate many $\hat x$ from many $y$ as long as the channel doesn't change at fixed complexity.
Complexity
You do the usual paper trick of ignoring the fact that you need to acquire $H$ first. That can, and quite possibly will, dominate the complexity of this (so don't be a bad scientist – figure out how you get your $H$ and how compelx that is, if you're writing some kind of publication).
So, assuming you already have the $H$, the operations you are doing per channel realization (i.e. only once every time $H$ changes, which probably doesn't happen very often):

calculating one half of the hermitian matrix $H^TH$ (no need to calculate the whole thing, it's redundant).
a QR decomposition of a hermitian matrix ($H^TH$)
calculating $\nu = Q^TH^T$

The operation you do per receive vector (i.e. for every symbol time step, typically, so pretty often):

Backsubstituting in $R\hat x$ to finally get $\hat x$.

The complexities per channel realization are:

$\frac{N^3}2$
depends very much on your QR decomposition algorithm (which depends on matrix size, and need for numerical stability, and your platform), but typically no worse than $N^4$
$N^3$

So that's a total of $C_H = \mathcal O \left(\frac 32 N^3 + N^4\right)$ multiply-accumulates.
The complexity  per receive vector:

$\frac12 N^2$

So, $C_y = \mathcal O \left(\frac12 N^2\right)$.
If we introduce a coherence time (measured in symbol durations) $T_C$ for which $H$ doesn't change, you get an overall complexity of
\begin{align}
C_{\text{overall}} &= \mathcal O \left(\frac {\frac 32 N^3 + N^4}{N T_C} + \frac1{2N} N^2\right)\\
 &= \mathcal O \left(\frac {\frac 32 N^2 + N^3}{T_C} + \frac12 N\right)
\end{align}
per receive symbol (i.e. per element of $\hat x$).
Bringing the constellation order into this
As said above, the constellation order only makes a difference to the complexity of the decider. However, we know under ML rule that the likeliest transmit symbol vector is the one geometrically closest to the likeliest transmit symbol vector – and since all elements of that are to be assumed independent from each other, that minimum can be found individually. The complexity of the decision thus only depends on the (SISO) constellation decoder you pick; nothing has been said about that, and there's immensely different implementations of these, so I will not allow myself to give you a complexity estimate here – I'll just allow myself that I can't see the formula as given by the paper appearing anywhere in my head.
It is a very valid question to ask whether it makes sense to apply the ML detection principle to individual elements of the transmit / receive vector. Unless you're in the very massive MIMO case, it makes sense to assume that there's channel coding on these symbols, and that they are, in fact, not stochastically independent from each other. So, classically, a MLSE-based (typically: Viterbi) soft-input channel decoder would be what you use, or some kind of message-passing based channel decoder.


Multiscale object detection generally consists in smoothing the image with larger and larger kernels, to check  which objects persist across different scales. So, zoom-out, mostly. You now can have access to the scale-space representation of the image, as shown below.

One  reason is that, without more knowledge, one does not zoom in beyond the pixel width. As illustrated (with humour) in the following picture.

Meanwhile, in certain conditions, subpixel resolution or superresolution is possible.


reduce harmonics produced by microcontroller

No, you need a low pass filter aka reconstruction filter.

In the general case your filter needs to reject everything at or above Nyquist = Fs/2 = 22.05 kHz. Typically it might have have a roll-off from say -3 dB at 20kHz to -80 dB at 22 kHz. This would be quite a demanding filter to implement with analogue components however, so if you are really only interested in generating a single tone at 18 kHz then you can use a less aggressive LPF which has a -3 dB point at say 20 kHz and a stop-band which begins at 35 kHz.


MFCC delta: How to compute delta, in real time, with no future frames?

As long as you train and evaluate the same features then you will be perfectly fine with this delay.

In practice first two frames of the feature vector for both $\Delta$ and $\Delta\Delta$'s are zeros. So in fact the coefficients are delayed by two frames. These can be calculated according to the formulas:

$$\Delta[t] = c[t]-c[t-2]$$
$$\Delta\Delta[t] = c[t-2]-2c[t-1] + c[t]$$

Notice that the scaling factors are omitted intentionally. It's only a linear scaling of the feature space and should not change the results. For example in the case of GMM's the Mahalanobis distance takes care of differences in variance across dimensions.


On the simplification using trigonometric functions

If you use the Euler's formula, you can simplify like this:
$$
[d]_{k,n} = \frac{\sqrt{2}}{N}\left( \cos{\left[ \frac{(k-1)(2n-1)\pi}{2N} \right]} e^{j\frac{2 \pi nk}{N}} \right)
$$

I think we can't simplify more.

PS: If you use your expressions of $A$ and $B$ and again the Euler's formula, you will get the same result.


Filters applied separately to right and left channel of stereo microphone yielding inaccurate results

I hope to distinguish between a voice input and musical input using a frequency cutoff.

The use case scenario is not exactly clear here. Are you trying to get the vacuum cleaner to be called to a location by responding to a command, or are you trying to get the vacuum cleaner to follow a person that is talking in general. (A specific person perhaps?).

The second use case scenario ("follow me as I am talking to you") is probably a "voice" recognition application and localisation could be added in two stages. In stage 1, a classifier that is trained to the voice of a specific person is fed frames from the sound card. Once the classifier detects that a frame contains a "voicing", then stage 2 is activated. In stage 2, you run a frequency analysis and try to locate the fundamental tone of the person's voice. Once you locate that, you then base your direction finding on the phase difference of that particular harmonic (or band). Obviously here, because the voice is processed in frames, the vacuum cleaner's reaction time depends on the length of the frame and their overlap. Very briefly, if you are processing 1 second long frames (with zero overlaps), then the robot makes 1 direction decision every 1 second.

The first use case scenario ("come towards me when I call you") is probably a pattern recognition problem. If the vacuum cleaner is to respond to a voice command, like "Robot, clean up", then, at the very least, you could use a matched filter approach. You can then rectify and integrate the output of the matched filter and pass it through a threshold. The combined output of this will be a pulse whose duration will be, approximately, equivalent to the duration of the voice command. Running this process on two microphones with some separation between them will give you two pulses (one on each channel) with a slight time difference from which you can then estimate direction.

A variation of this technique would be to call the robot via a whistle pattern (sounding a bit like R2D2 but obviously, at human speeds). The whistle pattern is attractive because:

You can filter your input to the expected range of the whistle (and it's very hard to whistle at 20kHz)
A whistle is very near to a sinusoid, therefore you can use very simple means of detection
Being a whistle and very close to a sinusoid, allows you to use the established techniques of phase difference to detect direction. In fact, for a given whistle pattern, you could choose to focus on the high frequency parts of it to get more accuracy out of the direction estimation.

So, given an audio stream which can even be processed on a sample-by-sample basis, establish a filter bank of filters that are atuned to a given frequency. Rectify and integrate their output (so that you essentially measure power at that band) and run that through a threshold. Once the user hits the right notes, you get a sequence of pulses and the times at which the outputs crossed the threhsold. Analysing the timings within one channel you can infer if the user just whistled the "secret pattern" and then analysing the timings across channels (or simply using a phase detector) you can infer direction.

Hope this helps.


How do you plot an SPL vs frequency curve for a given set of pressure data?

The steps will be, in this order:

apply window function (e.g. Hann or Hamming)
calculate FFT
calculate magnitude of FFT (sqrt(re*re+im*im))
convert magnitude to dB (20*log10(magnitude))
apply combined correction/calibration factors (window function compensation, FFT factor of N scaling if needed, factor of 2 for symmetric FFT, Pref calibration, etc) - this can just be a single dB correction value that you add to the dB magnitude values

If your Pref calibration data varies with frequency (e.g. microphone with non-flat frequency response) then it's a little more complicated than the above, but the same basic principles apply.

Note also that if you want weighted dB SPL (e.g. A weighting) then you'll also need to apply the relevant frequency dependent weighting corrections.


Can someone show the details of how to apply AIC for sinusoidal models to specific data?

The PDF (likelihood function) of the multi-tone estimation problem is:
$$
{\cal L}({\bf x}; {\bf a}, {\bf f},{\bf \phi}, {\sigma}^2) = \frac{1}{(2\pi \sigma^2)^{\frac{N}{2}}}
\exp\left( - \frac{1}{2\sigma^2} \sum_{n=0}^{N-1} ( x(n) - \sum_{p=1}^{P} {a}_p \cos(2\pi {f}_p + {\phi}_p) )^2   \right)
$$

The paper by Djuric quotes the AIC and MDL values as:

The code at the end attempts to:

Generate a 5-tone noisy signal.
Do the ML estimates of the amplitudes, frequencies and phases for a given value of $p$.
Find the (log) likelihood values for each value of $p$.
Calculate the AIC and MDL values from Djuric.

The plot for $p = 1 \ldots 60$ is:

Usually, the AIC over-estimates the model order. In this particular case, AIC gets 7 and MDL gets 5 (the correct value). The minima of each curve are plotted (as a o on the AIC curve and as a + on the MDL curve).

Djuric, P.M., "A model selection rule for sinusoids in white Gaussian noise," Signal Processing, IEEE Transactions on  (Volume:44 ,  Issue: 7 ), pp. 1744 - 1751.

scilab CODE ONLY BELOW

P = 5;
f = [ 0.1 0.12 0.231 0.333 0.478634];
phi = [0.329847 5.90324 3.09834 4.907983 1.984];
a = [ 1 1 2 2 3];
sigma = 1;

N = 100;

x = cos(2*%pi*[0:N-1]'*f + ones(N,1)* phi )*a';

xn = x + sigma*rand(N,1,'normal');

function [L,LL] = likelihood(data,ahat,fhat,phihat,sigmahat)
    N = length(data);
    datahat = cos(2*%pi*[0:N-1]'*fhat + ones(N,1)* phihat )*ahat';
    L = 1/(2*%pi*sigmahat*sigmahat)^(N/2) * exp(-1/(2*sigmahat*sigmahat)*sum((data - datahat).^2));
    LL = - log(2*%pi*sigmahat*sigmahat)*N/2 - -1/(2*sigmahat*sigmahat)*sum((data - datahat).^2);
endfunction

function [a,m] = aic(p,L,N)
    a = 3*p - log(L);
    m = 3*p/2*log(N) - log(L);
endfunction

function [ahat,fhat,phihat,regr] = sin_estimate(data)
    PAD = 10;
    N = length(data);
    DATA = fft([data; zeros((PAD-1)*N,1)],-1);
    [mx,ix] = max(abs(DATA(1:N*PAD/2,1)));
    ahat = mx/(N/2);
    fhat = (ix-1)/N/PAD;
    phihat = atan(imag(DATA(ix)),real(DATA(ix)));

    regr = data - cos(2*%pi*[0:N-1]'*fhat + ones(N,1)* phihat )*ahat
endfunction

pValues = [1:60];
likeValues = [];
loglikeValues = [];
aicValues = [];
mdlValues = [];
for p=pValues,
    signal = xn;
    a_est = zeros(1,p);
    f_est = zeros(1,p);
    phi_est = zeros(1,p);
    for pp=1:p
        [a_est(pp),f_est(pp),phi_est(pp),signal] =  sin_estimate(signal);
    end

    [l,ll] = likelihood(xn,a_est,f_est,phi_est,sigma);

    likeValues = [likeValues l];
    loglikeValues = [loglikeValues ll];

    [aicv,mdlv] = aic(p,likeValues(p));
    aicValues = [aicValues; aicv]
    mdlValues = [mdlValues; mdlv]
end

clf;
plot(aicValues);
plot(mdlValues,'r');
title('AIC and MDL plots');


Transform one-pole cascaded system into a biquadratic section

A single one-pole filter satisfies the following difference equation:

$$y[n]=x[n]+b\left(y[n-1]-x[n]\right)=(1-b)x[n]+by[n-1]\tag{1}$$

which corresponds to this transfer function:

$$H(z)=\frac{1-b}{1-bz^{-1}}\tag{2}$$

If you cascade two of these systems, the transfer functions are multiplied:

$$H^2(z)=\frac{(1-b)^2}{(1-bz^{-1})^2}=\frac{(1-b)^2}{1-2bz^{-1}+b^2z^{-2}}\tag{3}$$

The transfer function in $(3)$ corresponds to the following difference equation for the cascaded system:

$$y[n]=(1-b)^2x[n]+2by[n-1]-b^2y[n-2]\tag{4}$$

Consequently, implementing $(4)$ is equivalent to implementing $(1)$ twice.


Modelling Unwanted Signal in a LMS Adaptive Filter

The LMS and many of the variants of Adaptive Filters (In the Linear System context) work in the following settings (Intuitive):

You have access to 2 signals.
One signal is the result of the other one when a Linear System is applied.

This sounds really limiting, yet in practice it is powerful and flexible.

In the settings you mentioned the most known and similar problem is the Echo Cancellation model.

Pay attention that your model can't be formed as a Linear System.
Namely the connection between the clean signal and the corrupted signal can't be described using Linear System.


Fourier transform and anti-trasform--identity missing

Different ways of showing it, depends where you start. Are you willing to accept that the Fourier transform of $\delta(t)$ is $1$ and vice versa (i.e., $\int_{-\infty}^\infty {\color{red}1} \cdot {\rm e}^{\jmath 2\pi f t} = \delta(t)$)? If so, it's easy:

$$\begin{align}
r(\tau) & = \int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f \tau} S(f) {\rm d}f \\
& = \int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f \tau} \int_{-\infty}^\infty {\rm e}^{-\jmath 2\pi f t} r(t) {\rm d}t {\rm d}f \\
& = \int_{-\infty}^\infty \int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f (\tau-t)} r(t) {\rm d}t {\rm d}f \\
& = \int_{-\infty}^\infty r(t) \int_{-\infty}^\infty {\color{red}1} \cdot {\rm e}^{\jmath 2\pi f (\tau-t)} {\rm d}f {\rm d} t \\
& = \int_{-\infty}^\infty r(t) \delta(\tau-t) {\rm d} t \\
& = \int_{-\infty}^\infty r(\tau) \delta(\tau-t) {\rm d} t \\
& = r(\tau) \int_{-\infty}^\infty \delta(\tau-t) {\rm d} t \\
& = r(\tau)
\end{align}$$

Step 1: Insert (2) into (1). Note that the inner integration variable is a new one, different from $\tau$. I call it $t$.
Step 2: Pulling the first exp inside the integral, combining the exps.
Step 3: Changing integration order (PSD and ACF are absolutely integrable), pulling out what does not depend on the inner integration variable $f$.
Step 4: Using the fact that the inverse Fourier of a constant is a delta (think $\tau-t$ as one variable here, then it's the inverse Fourier of 1).
Step 5: Using the sifting property of the delta.
Step 6: Moving out the constant term
Step 7: Area under the delta is one.

Of course, step 4 is the critical one. If you don't buy it, you need a different, more fundamental/mathematical approach. This is more the engineering point of view I'm presenting here.

Regarding your reply with $\tau$ vs. $t$: What you say is not true. See, we're computing $r(\tau)$ via $\int_{-\infty}^\infty {\rm e}^{\jmath 2\pi f\tau} S(f) {\rm d}f$, which means it's an integral over frequency and the integration kernel depends on $\tau$. The function $S(f)$ does not depend on $\tau$ of course. Now, you are replacing $S(f)$ by the inverse Fourier transform of the autocorrelation, which is an integral over time. But it's a different time variable, which I called $t$. It must be, since it if were $\tau$ it would mean that $S(f)$ somehow depends on $\tau$.

Another way to think about it: The variable $\tau$ is the independent one. Don't forget that all our integrals are definite ones (from $-\infty$ to $\infty$). We sometimes drop that for laziness, but I added them now to be more clear. Now, the integration variables on the right-hand side are the ones we integrate over, hence the result cannot depend on it. Imagine something like $x(\tau) = \int_{-\infty}^\infty \int_{-\infty}^\infty X(f) g(\tau) {\rm d}f {\rm d}\tau$. This does not make sense as the right-hand side integrates over $f$ and $\tau$ (the result is a number) whereas the right-hand side depends on $\tau$. This is why we need a new time variable.


where we use estimation and detection in communication system?

Don't worry too much about defining these terms too precisely, because they are used in many contexts with slightly different meanings.
In very general terms, "estimation" is the calculation of a signal parameter, for example the phase, the mean, the PSD, etc. In other words, you have a signal, possibly noisy or distorted, and you want to find something about it.
"Detection" is concerned with whether you actually have a signal or not; think about radar, where you send a pulse, and wait to see if you detect the reflection. Since recovering the information from a PAM signal is based on the same principle (i.e. matched filtering), in digital communications it is common to say that transmitted symbols are detected.


Is Kalman Filter for LTV systems applicable on Hybrid systems?

You really have two questions -- one, is the Kalman filter applicable to linear hybrid systems, and the other is why isn't it used.
I can't answer your second question.
As for your first question (is it applicable) -- I think it could be made to be applicable.  In fact, I strongly suspect that there's some fairly obvious combination of a Kalman-Bucy filter for the continuous-time part of the system and a straight Kalman filter for the discrete-time part.
However
I've found that for a lot of real-world problems involving a discrete-time controller and a continuous-time plant, it's easier to just model the continuous-time plant as a discrete-time system as seen by the controller, optimize the control in discrete time, and then just sweep the fact that the plant is actually continuous time under the rug.
This is because if you're going to make a discrete time controller that actually works in practice, the most reliable way of achieving this end is to make sure that the frequencies of any significant plant dynamics are exceeded by the controller sampling rate by a healthy margin.  Do this, and any inaccuracies that you introduce by assuming that the plant responds step-wise to controller commands is less than any inaccuracies from the difference between your best-ever attempt at a plant model and what a real plant is going to do under real circumstances.
However (again)
If you're casting about for a PhD topic, there may be something in here, unless the approach was considered in 1952 and discarded for some good reason.


Calculate Q factor of a Low Shelf and High Shelf filter

it looks like this eq10q uses cascaded parametric EQs, initially set up to be bell, but switchable to other types (like shelving).

i wouldn't doubt that they're using the audio EQ cookbook.  if you're using that, there is a parameter called "shelf slope" or $S$ that is often set to 1, because that gets you the steepest slope without dips or lips or bumps.  so it's a straight incline or decline to the shelf level.  if the control is $Q$, instead, there is a relationship between $S$ and $Q$:

$$ \frac{1}{Q} = \sqrt{ \left( A + \frac{1}{A} \right) \left( \frac{1}{S} - 1 \right) + 2} $$

where

$$ A = 10^{dB_\text{gain}/40} $$

and $ dB_\text{gain} $ is the shelf gain in dB.

i would just plug in $S=1$ and get a $Q=\frac{1}{\sqrt{2}}$.


Most accurate method of Voice Activity Detection?

If you have more than 1 detectors, you may try blind separation based on the power spectrum difference of voice and background noise


How to deconvolve dependent part of signal from independent part?

First of all note that there is a certain ambiguity in the problem formulation because the approximation of $x(t)$ has to be divided into the part represented by $x^{\prime}(t)$ and the part represented by $f(y(t))$, and this division is not unique. The approach I chose is to start with a parametrization of $f(y(t))$, and approximate the given $x(t)$ as well as possible by this function. The remaining approximation error is then taken care of by $x^{\prime}(t)$. I use a least squares criterion for approximating $x(t)$ by $f(y(t))$.

I use the following parametrization for $f(y(t))$:

$$f(y(t))=a_1y^2(t)+a_2y(t)+a_3\tag{1}$$

with unknown parameters $a_i$, which need to be determined given the known functions $x(t)$ and $y(t)$.
Of course other choices for $f(t)$ are possible, but let's for illustration purposes use the above parametrization.

First you need to define a grid of values $t_i$ and set up an overdetermined system of equations:

$$x(t_i)=x^{\prime}(t_i)+a_1y^2(t_i)+a_2y(t_i)+a_3,\quad i=0,1,\ldots N-1\tag{2}$$

with some $N$ (which must of course be much greater than 3, the number of unknowns). Define now a $N\times 3$ matrix $\mathbf{Y}$ with rows $Y_i=[y^2(t_i), y(t_i),1]$, a $N\times 1$ vector $\mathbf{x}$ with elements $\mathbf{x}_i=x(t_i)$, and an unknown $3\times 1$ vector $\mathbf{a}$ with elements $a_i$. The optimal least squares solution for $\mathbf{a}$ is given by

$$\mathbf{a}=(\mathbf{Y}^T\mathbf{Y})^{-1}\mathbf{Y}^T\mathbf{x}\tag{3}$$

which in Matlab or Octave is simply implemented by the command

a = Y\x;

The values $x^{\prime}(t_i)$ are then simply given by

$$x^{\prime}(t_i)=x(t_i)-[a_1y^2(t)+a_2y(t_i)+a_3]\tag{4}$$

with the parameters $a_i$ as determined from (3). Note that $x^{\prime}(t)$ as defined by (4) is zero-mean. This is no restriction of generality because the parameter $a_3$ defining the constant part of $f(y(t))$ can as well be added to $x^{\prime}(t)$. This is part of the ambiguity inherent to the problem.

Using your example, the corresponding Matlab/Octave code would be something like

t = linspace(0, 100, 100);
y = sin(pi * t / 5) + t/50 + 1;
x = sin(pi * t/100) + 0.05 * y .* y;
x=x(:); y=y(:);
Y = [y.*y,y,ones(100,1)];
a = Y\x;
f = a(1)*y.*y+a(2)*y+a(3)*ones(100,1);
xprime = x - f;
plot(t,x,t,xprime,t,f)

From which you get the following plot ('x' blue, 'xprime' green, 'f' red):

As you can see, 'xprime' is zero mean, but his could easily be changed by adding the constant term determined by $a_3$ to 'xprime' (and subtracting it from 'f').

UPDATE:

So far no constraints have been imposed on $x^{\prime}$. If smoothness of $x^{\prime}$ is desired, this can be accomplished by adding more equations to the overdetermined system (2). These equations can be of the form

$$w\left[a_1(y^2(t_i)-y^2(t_{i-1}))+a_2(y(t_i)-y(t_{i-1}))\right]=w[x(t_i)-x(t_{i-1})],\\\quad i=1,\ldots,N-1\tag{5}$$

with some positive weight factor $w$ determining the relative importance of these additional smoothness constraints. If the matrix $\mathbf{Y}$ and the vector $\mathbf{x}$ are augmented by the values given in (5), the overdetermined system can be solved as before by solving equation (3). I computed the solution for $w=10$. Furthermore I added the constant term $a_3$ to $x^{\prime}$, which is an arbitrary choice (you might as well attribute the constant term to $f(y(t))$). The result is shown below (same functions and colors as in the above plot):

The plot was created in Octave by the following code:

t = linspace(0, 100, 100);
y = sin(pi * t / 5) + t/50 + 1;
x = sin(pi * t/100) + 0.05 * y .* y;
x=x(:); y=y(:);
Y = [y.*y,y,ones(100,1)];
w = 10;
y2 = y(2:100); y1 = y(1:99);
x2 = x(2:100) - x(1:99);
Y2 = w * [ y2.*y2 - y1.*y1, y2 - y1, zeros(99,1)];
Y = [Y;Y2];
X = [x;w*x2];
a = Y\X
f = a(1)*y.*y+a(2)*y;
xprime = x - f;
plot(t,x,t,xprime,t,f+a(3))


Reconstructed output mismatch for LTI system

The problem is that the formula $$H(f) = \frac{Y(f)}{X(f)},$$ where $H(f)$ is the frequency response, and $X(f)$, $Y(f)$ are the input and output, is valid only for frquencies $f$ where $|X(f)| \neq 0$.

Furthermore, if $|X(f_0)| \approx 0$, you're going to run into numerical errors, since $|H(f_0)| \rightarrow \infty$.

Your input signal has a bunch of places where it's zero, or close to zero. I wrote a little function that takes care of those cases:

function h = impres(x,y)

X = fft(x);
Y = fft(y);

for i = 1:length(X)
    if abs(X(i)) > 10
        H(i) = Y(i)/X(i);
    else
        H(i) = 0;
    end
end

h = ifft(H);

end

While not perfect, and could be tweaked, I get:

where the red line is the estimated output.

One more thing to be aware of: the equivalent in time of the IDFT is the circular convolution, not the linear convolution. In this case there is not too much difference, but you want to be as correct as you can in your code.


Is `xcov()`More Accurate Than `xcorr()` to Compute the Cross Correlation of Signals with Non Zero Mean Values?

When trying to measure similarity between signals we're basically building a metric.
When doing so we need to ask what we want to be sensitive about.
For instance, if you don't remove the DC Component (The Mean) and use something like an integral to measure similarity (Correlation / Convolution based) then you are sensitive the added DC component.
For instance, take 2 random white noises and add to each a DC component of 1e6.
Measure their similarity by correlation and compare it to 2 Sine Waves with amplitude of 1e-3 and a relative shift phase of 1e-6 [Rad].
In your eyes, which should be more similar? Which is indeed?
I'd say that in most cases in signal processing, similarity should be measured ignoring the DC component.


Eventhough it's not fully clear what sort of an environment you are referring to by the measurement set-up, in the most broadest sense of the term you need either a data acquisition interface or a sound card for the cheapest semi-alternative.

A data acquisition interface to a PC is what you mainly need. There are a bunch of alternatives. You must set your requirements (data rate, #channels, bit depths, ADC linearity, SNR etc.) to get the optimal match for the price. Note that the software interface to the acquisition hardware is another concern that makes big functional differences in the end.

The cheapest alternative to a professional data acquisition hardware can be a PC audio soundcard. They will record AC coupled and distorted signals. So no DC is captured unles you modify (if possible) by shorting the soundcard's ADC chip's input-DC blocking caps. There are a number of soundcard capturing software which also simulate an oscilloscope as well. Most typical PC (laptop) sound cards do provide 16 bits, 2 ch, 96 kHz sampling these days.

A logic analyser is almost useless from a DSP point fo view. Logic analysers are used for digital electronics such as logic gates, microprocessors, interface protocols etc. to analyse their bits, high and low values, ones and zeros.

An oscilloscope is again useful mainly for analog electronics diagnosis purposes. It's an indispensible hardware tool for any hardware project, but not exclusively required for data capturing, unless you want to observe in real-time what you are capturing.


The text:

  The measurements are influenced by gaussian noise with center in 0V and standard deviation equal to 3 Volts.

simply means that the Gaussian noise has $\mu = 0$ (the mean is zero) and $\sigma = 3$ (the standard deviation is 3 volts).

  Determine the model for the system, as well as initial values to implement a Kalman Filter

Without more information, a Kalman filter seems like overkill.

As @MBaz says, one model you could use is:
$$
r = V + n
$$
which has no time-dependent component. Making a Kalman-amenable model out of this might look like:
$$
V_{k+1} = V_k + m_k\tag{1}
$$
$$
r_k = V_k + n_k\tag{2}
$$
and we'd use the 14V nominal voltage to say $V_0 = 14$.

However, we are given no information about $m_k$ - how the true voltage might vary. We might also want to augment the state (perhaps the variation if $V_k$ depends on the derivative of $V+k$?) and so include a state transition matrix in front of $V_k$ in (1). If we did augment the state, we'd also probably want to add an output matrix on (2).

  Could you give me an insight about the initial covariance? Can it be reached from the information given?

Generally, I choose the initial state covariance to be VERY large. This is because I am generally given zero information about what it should be, as in this case. In this case, "VERY" might mean 30 (ten standard deviations).


Why we have to multiply the images with $(-1)^{x+y}$ to center the transform in frequency domain filtering?

Below 1D argumentation also explains the 2D case.

First consider the DTFT property for the pair $x[n] \longleftrightarrow X(e^{j\omega})$

$$ e^{j\omega_0 n} \cdot x[n] \longleftrightarrow X(e^{j(\omega - \omega_0)}) $$

Then recognise that $(-1)^n = e^{j \pi n} $ which yields:

$$ e^{j\pi n} \cdot x[n] \longleftrightarrow X(e^{j(\omega - \pi)}) $$

The effect is such that the spectrum $X(e^{j \omega})$ is centralized (shifted) into the $\omega = \pi$ frequency.

In 2D, the DTFT of the image is shifted into the central zone.


Convolution and Modulation Theorem with Octave Example

You are making the standard mistake of computing a linear convolution by performing a circular convolution which is implemented with an inverse DFT of wrong size.

Now assume that your signals $x[n]$ and $h[n]$ are of length $L_x$ and $L_h$ samples each. Then the discrete-convolution of $x[n]$ with $h[n]$ would produce the new sequence $$y[n] =  h[n] \star x[n]$$ whose length is $$L_y = L_x + L_h - 1.$$ This convolution is also referred to as the linear convolution of $x[n]$ and $h[n]$.

According to the well known theorem of DTFT the linear convolution of $x[n]$ and $h[n]$ is related to the multiplication of their DTFTs as: $$ h[n] \star x[n] \longleftrightarrow X(e^{j \omega})H(e^{j \omega})$$ where $H(e^{j \omega})$ and $X(e^{j \omega})$ are the theoretical DTFTs of $x[n]$ and $h[n]$.

The practical verification of this theorem is based on the multiplication of DFTs $X[k]$ and $H[k]$, however, which computes a circular convolution instead of a linear one. So one has to be careful in selecting the correct length for both forward and inverse DFT sizes.

In your implementation you have selected a shorter than necessary length so that there will be aliasing. Apart from specific instances aliasing should be avoid. The following MATLAb/Octave code shows the corrected version.

Lx = 12;           % Length of x[n]
Lh = 15;           % Length of h[n]
Ly = Lx + Lh -1;   % Length of y[n] = x[n]*h[n]

x = randn(1,Lx);   % sequence x[n]
h = randn(1,Lh);   % sequence h[n]

y1 = conv( x, h ); % y[n] computed in time domain by direct linear convolution
y2 = real( ifft( fft(x,Ly).*fft(h,Ly) , Ly) ); % y[n] computed by DFT method


Does a raised-cosine filter have infinite support?

An ideal raised cosine impulse is band-limited with cut-off frequency $f_c=(1+\beta)/(2T)$, where $\beta$ ($0\le\beta\le 1$) is the roll-off factor. Since it is band-limited, it must have infinite support in the time domain. The important difference with a sinc impulse is that its envelope decays much faster (for $\beta>0$, of course).

Note that even though an ideal raised cosine filter has infinite support, in a practical implementation you would use a filter with finite support (implemented as an FIR filter) that approximates the ideal raised cosine impulse in some interval around its maximum.


Maths behind why removing values from DFT is not a brick wall filter

A true brickwall filter would require an infinite impulse response in the time domain, since the inverse Fourier Transform of a rectangular window (brickwall frequency response) is a Sinc function. Since for the DFT the number of samples in the frequency domain is equal to the number of samples in the time domain, this would require you to have an infinite number of samples in the frequency domain and then require you to null an infinite number of samples.

For anything less, the impulse response in time becomes an aliased Sinc function which no longer corresponds to a true brickwall frequency response.

If you restrict your frequency space to just the DFT samples however, you will of course achieve a brickwall response following the rules that all samples outside of your rectangular window (in frequency) are zero. However in the implementation of a filter we are often concerned with the frequency response of signals that exist on a continuous frequency domain, in which case as described above you will not be able to achieve a perfect brickwall frequency response.


Image classification using SIFT features and SVM

If you could implement an SVM, you can quantize the features. :)

Typically the features are quantized using k-means clustering.  First, you decide what your "vocabulary size" should be (say 200 "visual words"), and then you run k-means clustering for that number of clusters (200).  The SIFT descriptors are vectors of 128 elements, i. e. points in 128-dimensional space. So you can try to cluster them, like any other points. You extract SIFT descriptors from a large number of images, similar to those you wish classify using bag-of-features. (Ideally this should be a separate set of images, but in practice people often just get features from their training image set.) Then you run k-means clustering on this large set of SIFT descriptors to partition it into 200 (or whatever) clusters, i. e. to assign each descriptor to a cluster. k-means will give you 200 cluster centers, which you can use to assign any other SIFT descriptor to a particular cluster.

Then you take each SIFT descriptor in your image, and decide which of the 200 clusters it belongs to, by finding the center of the cluster closest to it.  Then you simply count how many features from each cluster you have.  Thus, for any image with any number of SIFT features you have a histogram of 200 bins.  That is your feature vector which you give to the SVM.  (Note, the term features is grossly overloaded).

As I recall, there was a lot of work done concerning how these histograms should be normalized.  I might be wrong, but I seem to recall a paper that claimed that a binary feature vector (i. e. 1 if at least 1 feature from this cluster is present, and 0 otherwise) worked better than a histogram.  You would have to check the literature for details, and the details are important.

Edit: The Computer Vision System Toolbox for MATLAB now provides the bag of features functionality.


Sensor rise time compared to experimental step input rise time

So, when you remember how the step response, impulse response and frequency response are related, you'll notice that if you, instead of an actual impulse integral (i.e. a step) use something that is wider, then the frequency response simply gets windowed to lower frequencies.

In other words, and to put it as an information-gathering problem:

if the bandwidth of your excitation signal is limited (e.g. by a slow rise time of your step), then the information you can get about a system reacting to that is limited to the frequencies in the excitation.

It's often physically hard to make a large jump excitation. But maybe smaller jumps work, and you'd be better of

Generating a pseudorandom (band-limited) white signal (and save it!)
exciting your system with that
correlating the output of your system with the known white signal

The idea is this: you get the frequency response as the transform of the impulse response.

A white signal $s$ has the property ($*$ representing convolution here, ${}^*$ meaning complex conjugation, which is identity for real-valued signals)

$$s(t) * s^*(-t) = \delta\text,$$

i.e. you can "generate" an impulse $\delta$ by correlating a white signal with itself.

Knowing that your system has a step response and impulse response $h$, and knowing that convolution is commutative, we can do the following to find that:

$$\delta * h = s(t)*s^*(-t) *h = s(t)*h*s^*(-t)$$

In other words: instead of exciting your system with an impulse, which can be hard to do, to get the impulse response, we excite it with a white signal, which we can (withing the bounds of our time and bandwidth restrictions) produce digitally.

Then, we just correlate the output of the system with the same white signal, and get a result that is the same as had we put on an impulse of the same bandwidth as the white signal – but the white signal can be made arbitrarily long, so that very small amplitudes of excitation suffice to deliver the same energy as a really high, short pulse.

This is something that is commonly done to identify systems (especially: communication channels) in digital communications, but it applies to any LTI system.


Fourier series expansion of $\exp[-j2kA\sin(\omega t)]$

They're talking about the Fourier series expansion of the periodic function

$$f(\tau)=e^{-jx\sin(\tau)}=\sum_{n=-\infty}^{\infty}c_ne^{jn\tau}\tag{1}$$

which sometimes appears in the literature, e.g., when computing the spectrum of a frequency modulated sine wave. The complex Fourier coefficients in $(1)$ are given by

$$c_n=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(\tau)e^{-jn\tau}d\tau=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{j(-n\tau-x\sin(\tau))}d\tau=J_{-n}(x)\tag{2}$$

where the last equality follows from the integral representation of the Bessel function of the first kind of order $n$.

Combining $(1)$ and $(2)$ we get

$$f(\tau)=\sum_{n=-\infty}^{\infty}J_{-n}(x)e^{jn\tau}\tag{3}$$

With $x=2k_0A$ and $\tau=\omega_Dt$ you obtain the expansion used in the paper you quoted.


Save computational burden by calculating the output of decimate filter at output rate?

Is it possible to write a code, in a for-loop style, that only calculates the decimated data, and outperforms the performance of built-in functions?

Yes.
(the way you asked this, here would be where the answer is over, but it's not what you meant, sooooo let's continue)
"for loop style" isn't very specific, and any algorithm that implements some form of convolution will iterate over elements, but sure:
For decimation by a factor of $M$:
Decompose your filter in $M$ subfilters. Round-robin the input samples into the individual subfilters. Every $M$ input samples, calculate the output as sum of the outputs of the subfilters.
There's many publications on polyphase filters. harris' multirate signal processing is the go-to reference; but there's really many other references, too, see for example this slide deck from the Binghamton University (I'm not even sure the graphics are from that slide deck's author, my apologies if I'm wrong, but I feel I've seen the graphics before).
I recently added a different answer on the same topic, with a bit more background on the polyphase "Magic" that allows the lower rate.

The built-in filter functions provided by Matlab and a lot of realtime DSP libraries are optimized to have best computation and storage performance, but they don't provide the capability of calculating the filter output only at output sampling rate.

That's not true; any real-time DSP framework I can think of has decimating / resampling filters that do exactly that. (example: GNU Radio's decimating FIR filter function; Liquid DSP's firdecim, csdr fir_decimate, the list goes on: this is really a core functionality of every real-time signal processing toolkit.)
I don't know how efficient matlab's implementation is, but upfirdn with an interpolation ratio of 1 and a decimation > 1 should do exactly that.


Baseband Signals have zero imaginary component

You're right, but you've misread the sentence:
Baseband Data Communication means that there's no $f_c$.
For example, your USB cable to your keyboard, that's baseband communication. It can only carrier real-valued signals, so you only get to choose one dimension, not two, of signal freedom.
This is opposed to passband systems, where you can modulate the baseband signal up to a carrier frequency, which means your negative (relative to $f_c$) frequency spectrum doesn't have to have the same content as the positive.


What is about the circular convolution in OFDM

Assuming you have $N$ symbols to transmit encoded in block $k$,
\begin{equation}
 s(k)
 =
 \begin{bmatrix}
  s_1(k) \\
  \vdots \\
  s_N(k)
 \end{bmatrix}
\end{equation}
Performing $N-$FFT at the transmitter, you now have
\begin{equation}
 x(k) = F^Hs(k)
\end{equation}
It is here that you insert a cyclic prefix of size $D > L - 1$, where $L$ is the number of taps in your FIR channel, i.e. at the receiver you have
\begin{equation}
 r(k)
 =H_{ISI}\begin{bmatrix}
  \overbrace{x_{N-D+1}(k)} \\
  \vdots \\
  \underbrace{x_N(k)}_{\text{CP}} \\
  x_1(k) \\
  \vdots \\
  x_N(k)
 \end{bmatrix}
 +
 H_{IBI}
\begin{bmatrix}
  x_{N-D+1}(k-1) \\
  \vdots \\
  x_N(k-1) \\
  x_1(k-1) \\
  \vdots \\
  x_N(k-1)
 \end{bmatrix}
\end{equation}
where $\lbrace c_k \rbrace_{k=0}^{L-1}$ are the channel taps. Note that the channel matrix is decomposed into two parts, ISI (intersymbol interference) and IBI (interblock interference)
\begin{equation}
 H = H_{ISI} + H_{IBI}
\end{equation}
i.e.
\begin{align}
 H_{ISI}
 &=
 \begin{bmatrix}
  c_0 & 0 & \ldots & 0\\
  \vdots & \ddots & \ldots & \vdots \\
  c_{L-1} & \ddots &  & \vdots \\
  \vdots & \ddots &  & \vdots \\
  0 & \ldots & c_{L-1} & c_0
 \end{bmatrix} \\
 H_{IBI}
 &=
 \begin{bmatrix}
  0 & \ldots & 0 & c_{L-1} & \ldots & \ldots & c_1 \\
  \vdots & \ddots & &  &  \ddots &  & \vdots \\
  \vdots & \ddots & &  &  \ddots &  & c_{L-1}\\
  0 & \ldots & &  &  \ldots &  & 0\\
  \vdots & \ddots & &  &  \ddots &  & \vdots \\
  0 & \ldots & &  &  \ldots &  & 0\\
 \end{bmatrix}
\end{align}
If we remove the first $D$ samples of $r(k)$, we have
\begin{equation}
\tilde{r}(k)
=
 \begin{bmatrix}
  r_{D+1}(k) \\
  \vdots \\
  r_{N+D}(k)
 \end{bmatrix}
 =
 \begin{bmatrix}
  0 & \ldots & 0 & c_{L-1} \ldots & c_0 & 0 & \ldots & 0\\
  \vdots &  &  &  \ddots &  & \ddots &  & \vdots \\
   0 & \ldots &  &  \ldots & 0 & c_{L-1} & \ldots & c_0\\
 \end{bmatrix}
\begin{bmatrix}
  x_{N-D+1}(k) \\
  \vdots \\
  x_N(k) \\
  x_1(k) \\
  \vdots \\
  x_N(k)
 \end{bmatrix}
\end{equation}
After a little bit of manipulation, you can see that
\begin{equation}
 \tilde{r}(k)
 =
 C F s(k)
\end{equation}
where $C$ is a circular matrix that contains
\begin{equation}
 \begin{bmatrix}
  c_0 & 0 & \ldots & c_{L-1} & \ldots c_1
 \end{bmatrix}
\end{equation}
in its first row. Now, applying a DFT, we get
\begin{equation}
 y(k)
 =
 F
 \tilde{r}(k)
 =
 F
 C F^H s(k)
\end{equation}
Since $C$ is circular, then $F
 C F^H$ is diagonalizable because any circular matrix is diagonal in a Fourier basis, hence you get
    \begin{equation}
  y(k) = D s(k)
 \end{equation}
    where $D$ is diagonal containing the eigenvalues of $C$.


Radar type application for narrow band only?

These constraints absolutely exist. There are the norm! We could only wish in our wildest dreams to use as wide as bandwidth as we like. There are many areas in a radar system that place limitations on how wide the bandwidth can be and we'll go over a few straight forward ones. Mainly we're talking about limitations due to the antenna and waveguide as well as signal processing requirements.
Antenna and Waveguide
With antennas, especially arrays, their geometries can be defined in terms of wavelengths. Given that wavelength, and therefore frequency, we can yield the antenna's directivity pattern. It's important to stress that when you typically look at an antenna pattern, it is for a single RF frequency only.
Even simple pulses have some bandwidth around the RF frequency but are usually not too large. This allows us to make the narrowband approximation that the antenna pattern is the same for frequencies near the RF frequency, which is how the pattern was defined.
With an exceptionally wideband waveform the frequencies transmitted are now farther apart and so the antenna pattern itself can change significantly intra-pulse. These affects are more extreme at lower RF frequencies where desired bandwidths are well within the order of the transmitted frequency.
These effects are undesirable because (including but not limited to):

The effective radiated power is no longer the same during the pulse.
Because the antenna pattern itself is changing during the pulse, it distorts the waveform itself and you incur SNR degradation as well as undesired phase behavior which may lead to other issues, such as errors in angle measurements.

The waveguide itself, as well as other receiver/exciter components, need to support the bandwidth of interest and thus will require the use of certain technologies and processes to support it. This is a nice way of saying that you're going to be spending a lot of money and making sacrifices in size, weight, and power (SWaP).
Signal Processing
Assuming that your RF hardware is taken care of, there are issues with data collection and front-end signal processing. High bandwidth waveforms require high sampling rates. Introducing higher sampling rates has the following overarching issues

High sampling rate analog-to-digital converters (ADCs) can be prohibitively expensive...for a good one.
Depending on the bandwidths and pulse widths desired, you may have to deal with a large number of samples. In critical applications such as a fire control radar your computing resource (usually FPGAs/CPUs) need to be able to take care of common tasks such as pulse compression, FFTs, memory copies, etc. within a very strict time schedule.

These are just a few examples. I hope it's clear that the narrowband restriction is many times the norm when designing a system. Thankfully, technology is advancing where we can get away with higher bandwidths easier than we could in the past.


5 images per person is really very few.  In this situation it makes sense to use as much data for the training set as you possibly can.  This means for each person randomly select 4 images for training, and one for testing.  You should use this procedure several times so that you have a number of training set/test set combinations.

If you have access to the Statistics Toolbox for Matlab, check out cvpartition function.

If not, then you should have a loop over all persons, and for each person use the randperm function to generate a randomly shuffled list of indices of the images.  Then you can take the image corresponding to the first index and add it to your test set, and put the rest into your training set. Or you can put 2 images into the test set, or however many you like.  Generally, the training set should be larger than the test set.

You should generate the training and the test set several times, each time calling rng with a different value, to reset the seed of the random number generator.  If you generate 10 training/test set pairs, then you evaluate your algorithm 10 times and get 10 accuracy values.  From them you can compute error bars, or you can compare to other algorithms using hypothesis testing.

As to how you actually classify the images, that is very much a separate question, and the answer depends on what's in the images.  What exactly is on that image you've posted?


The thing you should keep in the back of your mind is that R Kalman published his seminal paper in 1962 (or there about).  He was a Control Engineer, interested in real systems at a time prior to widely available high fidelity numerical multi-physics modeling. The kinematic models you contrast, are derived from
$$
F=ma =m \frac{dv}{dt}.
$$
and $m$ is often a constant point mass.
  Planes and rockets consume mass for propulsion so we can see early on that there will be neglected terms.  $F$ should also include aerodynamic drag and for laminar flow, drag is roughly proportional $v^2$ which is not even close to the Navier-Stokes  equations.  So even a simple drag model is nonlinear and for a real time online algorithm a linear state space model is all he really had.  The solution to these oversimplifications is noise and small step sizes and prediction feedback.

Noise can also be actual random phenomena.  So a lot of non trivial modeling errors are lumped together in
$$
\frac{d}{dt} \mathbf{x}(t)= \mathbf{A}\mathbf{x}(t) + \mathbf{B} \mathbf{u}(t)
$$
where $\mathbf{u}(t)$ can be control forces, random forces, or psuedo modeling errors that we treat as noise. We also have a measurement set of equations that can also include noise, real or conceptual.
You can stick with this form and use a Kalman-Bucy Filter or you can convert the continuous state space model to a difference equation which is typically used in what is a Kalman Filter. To do that you need to integrate the equation just given.  I'm not going to go in detail because there are books and websites that show this but one start with calculating the matrix exponnential
$$\mathbf{F}= \mathbf{e}^{\mathbf{A} \Delta t} $$.  Some of your examples calculate the matrix exponential and some don't. You can often obtain a symbolic solutution for the matrix exponential or you can solve it numerically.

Getting back to your original question is that these models are simplifications of simplifications.  There may be poor choices taken but there really isn't a single correct discrete time 3D state space model. some may have $x_{k+1}=x_{k} + \Delta t v_{k}$ and some might include $\frac{(\Delta t)^2}{2} a_k$.  Real objects aren't point masses.  Angular orientation matters.  A lot of what is missing is treated as noise over small time steps.  When the nonlinearity is more severe, we switch to the extended Kalman Filter.  In the future, the state models will be real time multi-physics numerical routines.

The best model is the one that works best with you data.


I know this is an old question, but maybe this is useful to someone.
In CMSIS documentation, in Biquad Cascade IIR Filters Using Direct Form I Structure, says that the used equation is:

$$y[n] = b_0 x[n] + b_1 x[n-1] + b_2  x[n-2] - a_1 y[n-1] - a_2 y[n-2]$$

You need to look closely to coefficients $a_1$ and $a_2$. You need to change the sign.

So in your case, if matlab gave you this result (SOS Matrix)

1,  0,  -1,  1, -1.9974,   0.99743,
1,  0,  -1,  1, -1.9995,  0.99947,
1,  0,  -1,  1,  -1.9937,  0.99371,
1,  0,  -1,  1, -1.9982,  0.99817,
1,  0,  -1,  1, -1.995,  0.99499,

then you should change it to:

static float32_t filter_coeffs[] = {

1,  0,  -1,  1.9974, -0.99743,
1,  0,  -1,  1.9995, -0.99947,
1,  0,  -1,  1.9937, -0.99371,
1,  0,  -1,  1.9982, -0.99817,
1,  0,  -1,  1.995,  -0.99499,
};


What are high frequencies and low frequencies in a signal?

The high frequencies contribute to the fast varying parts of the signal(the sharp transitions),
while the low frequencies contribute to the slow variations of the signal in the time domain.

You might want to take a look over here:
http://cns-alumni.bu.edu/~slehar/fourier/fourier.html

Also, if you have the time you might want to take a look at this online course:
http://academicearth.org/courses/the-fourier-transform-and-its-applications
You won't regret it.


Find state space model from transfer function

So, assuming that the transfer function is between the output $Y(s)$ and the input $U(s)$, namely
$$
\dfrac{Y(s)}{U(s)} = G(s)\qquad\qquad(1)
$$
multiplying (1) by $U(s)\cdot(s^2-2s+1))$ yields
$$
s^2 Y(s) -2s Y(s) + Y(s) = sU(s) + U(s)
$$
Going back in the time domain we obtain
$$
\ddot y -2\dot y + y = \dot u + u\qquad (2)
$$
Now, we look for a realization of the kind
\begin{align*}
\dot x &= Ax+Bu\qquad (3)\\
y &=Cx
\end{align*}
with $x=(x_1,x_2)$ and
\begin{align*}
A &= \begin{pmatrix}0 & 1\\a_1 & a_2\end{pmatrix} & B&=\begin{pmatrix}0\\1\end{pmatrix}, & C&=\begin{pmatrix}c_1&c_2\end{pmatrix}\end{align*}
The next step is to find the values of $(a_1,a_2,c_1,c_2)$ for which (3) has the same input-output behaviour of (2). From (3) we have
\begin{align*}
 y&=c_1x_1+c_2x_2\\
\dot y&= c_1\dot x_1 + c_2\dot x_2 = c_1x_2 + c_2a_1x_1 + c_2a_2 x_2 + c_2u\\
\ddot y&= c_1\ddot x_1 + c_2\ddot x_2 = c_1a_1x_1+c_1a_2x_2 +c_1u+ c_2a_1x_2+ c_2a_2a_1x_1+c_2a_2^2 x_2 + c_2a_2u + c_2\dot u
\end{align*}
substituting into (2) yields
\begin{align*}
0&=(c_1+c_1a_1+c_2a_2a_1-2c_2a_1)x_1 + (c_2+c_1a_2+c_2a_1+c_2a_2^2-2c_1-2c_2a_2)x_2\\
&+(c_1+a_2c_2-2c_2-1)u + (c_2-1)\dot u
\end{align*}
since that equality must hold for all $(x_1,x_2,u,\dot u)$ that' equivalent to ask
\begin{align*}
c_1+c_1a_1+c_2a_2a_1-2c_2a_1&=0\\c_2+c_1a_2+c_2a_1+c_2a_2^2-2c_1-2c_2a_2&=0\\
c_1+a_2c_2-2c_2-1&=0\\c_2-1&=0
\end{align*}
and a solution is
$$
(a_1,a_2,c_1,c_2)=(-1,2,1,1)
$$
Therefore
\begin{align*}
\begin{pmatrix}\dot x_1\\\dot x_2\end{pmatrix} &=\begin{pmatrix}0&1\\-1&2\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix} + \begin{pmatrix}0\\1\end{pmatrix}\\y&=\begin{pmatrix}1&1\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}
\end{align*}


Why is a negative exponent present in Fourier and Laplace transform?

Matt is correct that the sign is convention.  I think that there is a reason for it beyond that though.

If we look at complex frequencies in the complex plane, they look like a constant vectors that rotate in one direction or another.  Positive frequencies rotate counter-clockwise, negative frequencies rotate clockwise, and "0 Hz" frequencies don't rotate at all.

The Fourier transform has a negative sign to intentionally rotate in the opposite direction as the frequencies that they are "looking" for.

The reason for the opposite rotation is that when the two frequency vectors are multiplied, their phases will repeatedly cancel out, so when the results are summed together there will be a massive vector due to all of the individual vectors lining up.

$$
X(f) = \sum\limits_{n=0}^{N-1}x(n)e^{-j2\pi kn/N}
$$

This is how the Fourier transform "looks" for frequencies.  If the two frequencies are the same or "close" (how close they need to be depends on the length of the DFT) they will line up well and cause a massive response in the summation.  I have showed how this works for the discrete Fourier transform (DFT), but the exact same reasoning applies to the continuous transform.

Hopefully this explains why the Fourier transform would want the vectors rotating in the opposite direction.  To be perfectly honest I don't know the Laplace transform well enough to give solid reasoning for its negative sign.  Since the two transforms are very closely related though (the Laplace transform being a generalization of the Fourier transform), I assume that it is for similar reasons.


The range resolution is normally some fraction bandwidth of the pulse. The fraction will depend on how you measure it i.e. null to null, 3dB points, 6 dB points etc. Also the windowing (Hamming, Kaiser) used will reduce the resolution i.e. increase the number. For a chirp pulse with the bandwidth you stated the null to null width is given by

$$ \delta r = \frac{c}{2BW}=\frac{3\cdot 10^8}{2\cdot  40\cdot 10^9 }=3.75 \cdot 10^{-3} \textrm{ m } $$

You can reference the following links:

High Range Resolution Techniques
Understanding Millimetre Wave FMCW Radars,
Radar Basics


Why bother impluse invariance, bilinear transform?

Simply put, the digital implementations of the IIR filters (bilinear) are much cheaper in terms of computation time and memory.  An FIR or FFT takes much longer to process than an IIR, and requires much more memory.

An FIR requires you have a buffer of the length of the FIR, and another buffer to hold the FIR values.  For every sample, you must convolve these two buffers.  Since FIR filters are often very long (hundreds of taps) it can be difficult to get enough memory or processing power in a small embedded system.

Even in desktop systems where you have gigabytes of RAM and gighertz of processing speed you can run more filters using IIR than you can FIR filters.

On a desktop PC, your main limit would be processing speed rather than memory.  In an embedded system, you often have limits on memory and processing speed.

Using an IIR implemtation let you get away with cheaper hardware or to get more out of the hardware you have.

An IIR filter can be implemented with a handful of multiply and add steps for each sample.  An FIR requires as many multiply and add steps as the length of the filter.  Using the FFT to implement the FIR will save you some computation, but it will still take much more processing and ram than the IIR.

Another thing is that an IIR filter with a sharp cutoff isn't much more complicated than one with a shallower cutoff.  An FIR requires some minimum of taps to work at all (depending on sampling rate and cutoff frequency,) and making it sharp and clean requires a greate many more taps than that minimum.


If you want to generate a very short pulse, it is not usually a problem to have it at any time you wish to have it. That is because modern audio interfaces (sound cards) have an interpolating lowpass filter in the output path, usually with an impulse response that without latency looks like the sinc function:

Figure 1. Sinc function, shown truncated at the ends and critically sampled with the sample number shown in the horizontal axis.

If you play back audio data that is sampled from a time-shifted sinc function with the peak somewhere between sample points, the interpolating lowpass filter will reconstruct the function from the samples so that the peak is in the correct position. Some audio interfaces use a low-latency variant of the lowpass filter that does not have a symmetrical sinc impulse response. Disregarding latency, it has just as good inter-sample timing capabilities if you give it samples of a time-shifted sinc function.

If you need a more rectangular pulse, you could add an analog comparator to the output of the audio interface. You can change the pulse length by varying the voltage level you compare to.

First test if you get that sinc shape from your hardware, like MBaz suggested.


How do you calculate the group delay in samples of a Biquad Filter at a given frequency from the coefficients?

Pick two different frequencies slightly above and below your target frequency. The choice depends on what numerical precision you have available and what your sample rate is.
Calculate the z-transform at both frequencies
Calculate the phase at both points, make sure to check for "phase wrapping"
Subtract the phases and divide by the frequency difference
Scale properly and multiply with $-1$


The short answer is: You can implement a lowpass filter by sampling its prototype continuous time impulse response $h_c(t)$, but you cannot implement a highpass filter by sampling its prototype impulse response, due to highpass filter being  non-bandlimited i.e. there will be unavoidable aliasing due to the definition of the highpass filter. But read on.

The technique of obtaning a discrete time filter impulse response $h_d[n]$ from a continuous time prototype filter by sampling its impulse response $h_c(t)$ as given by $$h_d[n] = T_s h_c(nT_s)$$ is called as the method of impulse-invariance , where $T_s$ is the sampling period.

This operation of obtaning a discrete time filter from a continuous time prototype is analogous to the more traditional operation of signal sampling. In fact their implications are the same: There will be aliasing at the resulting discerete time spectrum of the converted filter if either 1-The continuous time filter is not bandlimited or 2-The sampling rate is not adequate, given by the Nyquist rate as applied to sampling of bandlimited signals.

A continous time lowpass filter by definition, is always bandlimited, and we can always, therefore, find a small enough $T_d$ which will avoid aliasing at the resulting discerete time spectrum of the filter being converted.

On the other hand, a continuous time highpass filter is going to a have an infinite bandwidth and as such the sampling period required is $T_s = 0$ , i.e. it cannot be sampled.

However, if you will be carefull enough, and if you can be sure about the bandwidth requirements of every stage in your signal processing chain, then may be you can replace a highpass filter by a bandpass one, when such is the case, then you can still apply the technique of impulse invariance to obtain a discrete time band pass filter from a continuous time bandpass filter prototype, by choosing a suitable sampling period $T_s$

Note that, the method of impulse invariance produces an exact linear mapping from continuous time to discrete time frequency spectrum. It will preserve all the characteristics of the continuous time filter (unlike the method of bilinear transform which warps the spectrum)


Where do you recommend getting "standard" images for testing algorithms?

Research communities don't "trust"; they take a set of data and consider its properties (sample size, representativity, quality, systematic biases...) and then pick that. There's no "universally trusted" thing, as a dataset that's perfect for one purpose might be very wrong for another.
So, you'll have to work this from the other side up: pick a set of images that well represent what you want to do! Is it improving the quality of compressed photographs of cattle? Different dataset then you'd use to assess a denoising algorithm for architectural photography!
If you want to do that works on most photographic content, for example, you'd probably end up using the same dataset as the JPEG used to develop their compression standards (of the same name); not only has that been proven to actually be a useful set, but also it allows a fair comparison to established methods that use the same.


Calculating similarity between two signal plots

Why not just use something like the relative "error" between the two?

For example, if your frequency magnitude responses are $G_1$ and $G_2$, then calculate:

$$
ERR = \sum \left | G_1(n) - G_2(n) \right|^2
$$

and then normalize with respect to the reference, $G_1$:

$$
NORMALIZED = \sum \left | G_1(n) - G_2(n) \right|^2 /  \sum \left | G_1(n) \right|^2
$$

Correlation is also a way to go, but it might show some cases where the same shape happens, but at very different frequencies... which might not be what you want.


How to distinguish signal and noise eigenvalues of autocorrelation matrix for MUSIC algorithm?

This problem is defined as estimating the dimension of a model. I would recommend use something like  Bayesian Information criteria (BIC) or Minimum description length (MDL).