Designing an FIR filter in Matlab with differing number of precursor taps and postcursor taps

There is no reasonable way to explicitly choose the number of taps before and after the peak. The reason is simple: the arbitrary magnitude response design results in a linear phase, and, consequently, the impulse response is symmetrical. If you specify a desired complex frequency response in terms of magnitude and phase, the location of the peak is an outcome of the design process. The smaller the average delay (implicitly specified by the desired phase response), the further to the left the peak will occur.

Take as an example the design of an FIR low pass filter with $21$ taps. The figure below shows $3$ designs. The difference is the desired phase response. The top figure shows the impulse response (filter taps) for a linear desired phase response, the middle one for a non-linear phase response with a lower pass band delay than the linear phase response, and the bottom figure shows the impulse response of the corresponding minimum phase filter. The location of the peak is entirely determined by the specified phase response.


What is the minimum phase shift for Quadrature Amplitude Modulation (QAM)?

The minimum phase shift in 16-QAM, 64-QAM, 256-QAM etc is 0 (same value regardless of whether you are using degrees or radians!) and occurs when the symbol transition is from one constellation point on some radius vector to another constellation point on the same radius vector.  If we take the constellation points as having coordinates $(i,j)$ where $(i,j) \in \{-3, -1, +1, +3\}$ in 16-QAM, then $(1,1)$ and $(3,3)$ lie on the same radius vector (line of slope 1 through the origin) and so the RF phase remains unchanged in the transition from one constellation point to the other. Obviously, the same example will also work for 64-QAM, 256-QAM, and 1024-QAM.

What is the minimum non-zero phase shift?

Well, for 16-QAM, the constellation points $(1,3)$ and $(3,3)$ are distance 2 apart, and together with the origin, form an obtuse triangle with sides of lengths $2$, $\sqrt{10}$, and $\sqrt{18}$. From the standard formula
$$c^2 = a^2 + b^2 - 2ab\cos C$$ relating the sides $a$, $b$, $c$ of a triangle opposite to angles $A$, $B$ and $C$, we get that
$$2^2 = 10 + 18 - 2\sqrt{180} \cos C$$
giving $\arccos \frac{2}{\sqrt{5}}$ as the minimum nonzero phase shift in 16-QAM.
Similar calculations for 64-QAM, 256-QAM and 1024-QAM as left as an exercise for the reader.


Ideal Discrete Time Hilbert Transformer with constant group delay

There's not much need for complicated calculations. If you have the zero-delay frequency response $H(e^{j\omega})$ as defined in the first equation of your question, the frequency response with a linear-phase (constant delay) is given by

$$H_d(e^{j\omega})=H(e^{j\omega})e^{-j\omega\tau}=e^{-j\left(\frac{\pi}{2}\textrm{sgn}(\omega)+\omega\tau\right)},\qquad \omega\in(-\pi,\pi)\tag{1}$$

where $\tau$ is the group delay in samples.


What is the difference between upsampling and up conversion (done during signal transmission)? Also is upconversion different from modulation?

Your first sentence is not correct.  The terminology of 'sample rate conversion' operations is often misused.  I suggest you replace your first sentence with, "Interpolation is when you increase the sample rate of a DISCRETE signal sequence."  I'd also replace your second sentence with, "Modulation is when you, by multiplying with a carrier frequency, shift a signal's spectral components up in frequency by a fixed amount." (The term "modulation" is typically used when talking about continuous signals.)  As far as I understand the terminology, "up-conversion" is when you shift the spectral components of a discrete sequence up in frequency by a fixed amount. [-Rick Lyons-]


Filter H(z) manipulation frequency response changes

It is just like ordinary functions manipulation.
You have to remember that z and the frequency are linked with $z = exp(j \omega \cdot T)$
$H(z^2)$ you are doing something equivalent to $H(2 \omega)$ on the continuous frequency domain, i.e. you are narrowing the filter band, and thus stretching its time response.
On the second example, you have to remember that $H_1(z) H_2(z)$ is the response of two cascading filter, that the product of the Z transforms is the convolution of the the time responses.
The filter you are suggesting $0.5 (1 - z^{-1})$ will be a sort of derivative estimate because a signal with response $x[n]$ will be transformed to $(x[n] - x[n-1])/2$ and this is a high-pass filter. You can also think of it as the inverse of the step response, you can spend some time contemplating a table of pairs of signals and their Z transform, for instance this and slowly but certainly you will build some intuition.


ROC of Z Transform of $x(n) = 2(3)^nu(-n)$

Using the definition, and letting $k=-n$:

$X(z)=2 \sum_{n=-\infty}^03^nz^{-n}=2\sum_{k=0}^{\infty}\big(\frac{1}{3}z\big)^k$

Now we can use the geometric series (https://en.wikipedia.org/wiki/Geometric_series) formula only if we require $|\frac{1}{3}z|<1$, or equivalently $|z|<3$


how to calculate the dc component of an energy signal?

An energy signal $f(t)$ satisfies

$$\int_{-\infty}^{\infty}|f(t)|^2dt<\infty\tag{1}$$

i.e. it has finite energy. If $f(t)$ had a non-zero DC component $f_0$ defined by

$$f_0=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{-T/2}^{T/2}f(t)dt\tag{2}$$

the integral in (1) wouldn't converge. So an energy signal must have a zero DC component.

Note that there is some confusion about the term "DC component". Some people would call $F(0)$, i.e. the value of the Fourier transform of $f(t)$ at frequency $\omega=0$, the "DC component" of $f(t)$. However, note that $F(0)$ given by

$$F(0)=\int_{-\infty}^{\infty}f(t)dt\tag{3}$$

can be non-zero (but finite), even if (1) is satisfied. One example of such a signal satisfying (1) but with $F(0)\neq 0$ is a sinc pulse.

Also have a look a this answer to a related question.


The Hilbert-Huang signal decomposition possesses a non-adaptive part (the Hilbert transform) and an adaptive, or data-driven part: the decomposition of  a signal as a sum for Intrinsic Mode Functions (IMF), that are obtained by a deflation process, using extrema and envelopes of the signal, recursively.

Thus, IMFs are adaptive with respect to the signal: if you change the value of one sample, the IMFs can be drastically different. For some classes of processes, namely fractional Brownian motions, the
Empirical mode decompositions [can act as] as data-driven wavelet-like expansions. In other words, under some conditions, an adaptive algorithm may behave like a non-adaptive one.

There is one code for Marginal Hilbert Spectrum on Matlab Central.


Butterworth filter approximation: derivation and output poles

Note that the Butterworth filter is only one of many possible approximations of an ideal low pass filter. It is a Taylor series approximation (at $\omega=0$) of the desired response, and, consequently, the magnitude response of the Butterworth filter is maximally flat at $\omega=0$.

Approximating the squared magnitude response (instead of the complex frequency response, or the magnitude response) has the advantage that the function to be approximated is an analytic real-valued function of a real-valued variable, which simplifies the approximation problem.

We have

$$|H(j\omega)|^2=H(j\omega)H^*(j\omega)=H(j\omega)H(-j\omega)=H(s)H(-s)\Big |_{s=j\omega}\tag{1}$$

because $H^*(j\omega)=H(-j\omega)$ for real-valued systems, i.e., for systems with a real-valued impulse response, or, equivalently, with real-valued transfer function coefficients.

Since for a Butterworth low pass filter

$$|H(j\omega)|^2=\frac{1}{1+\omega^{2N}}\tag{2}$$

we have

$$H(s)H(-s)=\frac{1}{1+(s/j)^{2N}}=\frac{1}{1+(-s)^{2N}}\tag{3}$$

From the squared magnitude function we need to derive the transfer function $H(s)$, which is the function that we actually want to implement. From $(3)$ we know that the poles of $H(s)H(-s)$ (i.e., the zeros of the denominator of $(3)$) lie one a circle in the complex plane. Since we can only implement a causal and stable transfer function, we choose the $N$ poles in the left half-plane for $H(s)$, and consequently, the mirrored poles in the right half-plane are the poles of $H(-s)$.

Now we know $H(s)$ (i.e., its numerator and denominator coefficients), and, after having chosen a cut-off frequency and a filter structure, we can compute the values of the passive elements necessary to implement the given transfer function.


IIR oscillator produces stable output in MATLAB, unstable on chip, how do I tune it?

A first thing: I think there is a $\pi$ missing in your conversion between frequencies and angular frequency. Your impulse response shows a period of about 90, while your period should be about 29 samples - there's clearly something wrong here.

What you are seeing on your target hardware is due to the difference in precision between matlab's implicit use of double precision (64-bit), and the use of 32-bit for the routines you chose to use.

Even if this is not very apparent through the transfer function or the biquad structure, what this "filter" with no input does is simply rotate its initial conditions. With 32-bit precision, there is no guarantee that the magnitude of the input vector will be preserved by rotation - it can either grow (instability) or decay.

If you intend to use this as an oscillator (that is to say, if the input will always be null except for the initial conditions):

Regularly - say every 100th sample - mess with your state variables to readjust the norm of the vector to 1.
Or use another simpler and more stable technique like a phase accumulator and a sine lookup table.


what is the relation between a minimum phase system and a FIR system?

The only FIR system that is both minimum phase and linear phase is where the output is the same as the input: y[i] = x[i]


Say you want to transmit complex symbols $s_0,s_2,\ldots,s_{N-1}$ using frequency-division multiplexing. In the "bank of modulators" approach, you  would assign a carrier frequency $f_k$ to each symbol and create the signal
$$x(t) = \sum_{k=0}^{N-1} s_k \exp(j2\pi f_kt).$$

Contrast this with the definition of the inverse discrete Fourier transform:
$$x[n] = \sum_{k=0}^{N-1} s_i \exp(j2\pi nk/N).$$

So, if you interpret the symbols as the magnitude and phase of each carrier, and choose the carriers as $f_k=k/N$, then you'll see that you can generate $x(t)$ simply by converting $x[n]$ to an analog signal.

Note that:

I've explained how to transmit one block of symbols; extending this to multiple blocks has a slight complication and requires inserting a cyclic prefix to $x[n]$.
I'm assuming a rectangular pulse shape, which results in sinc-shaped carriers in frequency. The beauty of this method is that the carriers are orthogonal even if they overlap. In practice, this means that you can recover the symbols $s_k$ by performing the DFT of $x(t)$.
It is perfectly possible to use other pulse shapes with this method.


If you can reliably determine which samples are on a rising portion and which are on a falling portion, then perhaps you can do a linear regression on each (near)monotonic segment, and solve for the intersections of the computed regression lines.

If you know the non-linearity equation, then you can do a regression against that curve instead of using linear regression, and compute or numerically solve for the curve intersections, if needed.

This regression method might also help with error or robustness estimates, if needed.


Inverse Laplace transform of two-sided and one-sided Laplace transform

The inversion formula is the same for both types of transforms:

$$f(t)=\frac{1}{2\pi j}\int_{\alpha-j\infty}^{\alpha+j\infty}F(s)e^{st}ds\tag{1}$$

The difference is in the choice of the constant $\alpha$. The line $\textrm{Re}\{s\}=\alpha$ must be inside the region of convergence (ROC). For causal functions (i.e., functions for which $f(t)=0$ for $t<0$), the ROC is to the right of the pole with the most positive real part, whereas for non-causal functions, the ROC is a vertical strip between two poles.


Is it possible to apply inverse Fourier transform (I-FFT) to the following image?

Suppose, I have this image in my hand and nothing else. Clearly, this is a Magnitude-plot of some unknown image.

Sounds like you're losing the phase information in the signal; you only have the pixel values, correct?

  Is it possible to apply an Inverse Fast Fourier Transform (I-FFT) operation to recover the original Grayscale image from this image?

The important information in images tends to edges -- edges contain most detail in the image. This information is carried in the phase information -- when decomposed into sinusoidal signals, a square wave can have arbitrary location only if the waves contain phase information.

Additionally, real world images tend to look roughly the same in the frequency spectrum when only plotting magnitude. These images tend to have a large (in magnitude) low frequency components with (much) smaller high frequency components. It's reasonable to say that a significant portion of the entropy (or information) is carried in the phase information.

With a naïve inverse Fourier transform on the values obtained from the image, it is not possible (at least by experiment) to recover the original signal. However, the phase only inverse Fourier transform (where we divide by the magnitude to only preserve the phase).

The data behind the image was generated with

x = rgb2gray(astronaut())
X = fft2(x)
X_mag = np.abs(X)  # magnitude only
X_pha = X / (np.abs(X)) # phase only
y_mag = ifft2(X_mag)  # both y_mag and y_pha have max(y.imag) < 1e-15
y_pha = ifft2(X_pha)  # (some Fourier transform theory here + floating point precision)


Gaussian FIR filter with no multipliers?

This is an approximation, but you can make it as good as you like.

   Just use a cascade of several filters with rectangular impulse responses. In the simplest case this would be a two-tap filter. This works because of the central limit theorem. However, you will need to scale, because otherwise your resulting impulse response may become too large. The scaling could be done by bit-shifting.


Normalization factor in the convolution theorem

The choice of the normalization factor is just a matter of convention. Note that the specific correspondence between convolution in the time domain and multiplication in the frequency domain with a scaling of $\sqrt{2\pi}$, as shown in your question, applies only to the unitary definition of the Fourier transform with angular frequency as the independent variable in the frequency domain:
$$F(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(t)e^{-j\omega t}d\omega\tag{1}$$
If we use the non-unitary definition (as is common in the field of signal processing)
$$F(\omega)=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}d\omega\tag{2}$$
then the correspondence is simply
$$\mathcal{F}\{(f\star g)(t)\}=F(\omega)G(\omega)\tag{3}$$
The same is true if we use frequency instead of angular frequency.
For the discrete-time Fourier transform (DTFT), we have exactly the same choices.
For the discrete Fourier transform (DFT), where we deal with samples in both domains, the same holds for the normalization factor $N$, where $N$ is the DFT length. In the field of signal processing, the DFT is commonly defined as a non-unitary transform:
$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{4}$$
The inverse transform is then given by
$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{5}$$
The unitary DFT definition is
$$X[k]=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{6}$$
and the corresponding inverse transform is
$$x[n]=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{7}$$


Why do we need deterministic measurement matrices in compressed sensing?

As far as I know there are two reasons:

In sensing part: For practical implementation, usage of random matrices
is hard, so people try to come up with simpler matrices that are fixed,
this is thought to be helpful for designing efficient hardware.
In reconstruction part: To find deterministic matrices which their
structure can be exploited for faster signal recovery. Using a
structured matrix, might pave the way for efficient recovery. The
goal is to achieving something like FFT for compressive recovery.


How do I calculate peak amplitude of the signal components after zero padding and FFT?

It looks to me like your frequency resolution is sufficient and you are getting the correct value of the peak magnitude in your FFT but your scaling factor is incorrect.

You are right in dividing the magnitude of your FFT peak by N (1024 in your case). This number however, divides the power across the positive and negative sides of the frequency spectrum.

If you only want to look at one side of the FFT, you should multiply your peak magnitude by 2/N to get the correct amplitude that you are expecting. See this answer for another example.

e.g.

512 * 2 / 1024 = 1 (the answer that you were expecting)


The link in my comment suggests the following:

  I would recommend first extracting the envelope by using either a halfwave
  rectification (i.e., replace all the negative values in the time waveform
  with zeros) or a Hilbert Transform and then lowpass filtering the waveform
  at around 50 Hz (the lowpass is optional if you only care about the 4 Hz
  component).   Once you have the envelope, simply do an fft (making sure
  that your frequency resolution is at least 1 Hz) and look for the energy
  around 4 Hz.  If any of these steps don't seem obvious too you, let me know
  and I can send you some Matlab code.

So, what that means is you need to take your speech $s[n]$ and find the envelope:

$$
e[n] = \left | s[n] + j{\bf H}\left[s[n]\right] \right |
$$

where ${\bf H}$ is the Hilbert transform of your speech signal.

Then, take the FFT and look at the bin at 4Hz.


In order to accomplish a particular specified filtering task, the FIR filter order will almost certainly be much higher than the IIR filter order.  For instance a low-pass filter that requires a particular "sharpness" in cutoff, the transition region between the passband (at low frequencies) and the stopband (at high frequencies).  What a 4th-order IIR filter can do might correspond to what a 20th or 25th-order FIR can do.  The 4th-order IIR requires 10 coefficients and 4 states.  the 20th-order FIR requires 21 coefficients and 20 states.

Feedback can be useful.


Your reasoning is correct and we arrive at the same conclusion.

The system's output is given by

$$y[n]=x[-n+2]\tag{1}$$

Let

$$x_2[n]=x[n-k]\tag{2}$$

The corresponding output is

$$y_2[n]=x_2[-n+2]=x[-n-k+2]\tag{3}$$

From $(1)$ the shifted output is

$$y[n-k]=x[-(n-k)+2]=x[-n+k+2]\tag{4}$$

Since $(3)$ and $(4)$ are not identical, the system is time-varying.


For an input with limited amplitude, the maximum amplitude gain of a filter can be calculated as the absolute sum of the impulse response
Scaling with the inverse of the absolute sum, will guarantee that your max output amplitude will never exceed the max input amplitude
In practice, this does NOT work well: It's a very conservative scaling that often results in bad signal to noise ratio and the input signal that would actually trigger the max amplitude gain is extremely unlikely to occur. Scaling this way leaves a lot of headroom on the table.
A better approach is to look at your worst case sine wave gain. Sine waves (or similar signal) do happen in music and audio a lot, so that's a much more likely scenario. Especially in your case, where the max is around 200 Hz.
That still may not result in acceptable signal to noise ratio: in this case you can "eye ball" the gain at the most important frequencies and use a limiter or soft-clipper to manage the occasional overshoot.

In your particular case, I'd probably go with -10dB scaling with a limiter (provided, it's implemented in floating point). At the end of the day this is a trade off between managing clipping/overshoots, efficiency, signal to noise ratio and max output capability. It needs to be optimized within the specific requirements of your application


How do MATLAB and/or Python treat $2^n$ samples rule in FFT

Matlab doesn’t require the number of samples to be a power of 2. You can read more about it in Matlab FFT help page.
It will handle any length of signal and will not truncate or pad the length unless you specifically say so.
BTW in the help page it says that Matlab FFT is based on FFTW. In the FFTW documentation it says:

  FFTW is best at handling sizes of the form 2^a 3^b 5^c 7^d 11^e 13^f, where e+f is either 0 or 1, and the other exponents are arbitrary. Other sizes are computed by means of a slow, general-purpose algorithm (which nevertheless retains O(n log n) performance even for prime sizes).


signal uncorrelated with delta in the origin

It has a delta at the origin because you are multiplying your time series with itself (i.e. there is no time delay).  Any non-zero time series correlated with itself will give a positive value.


Plotting the spectrum of a continuous sinusoid function

Given your signal $x(t) = 5\cos(800\pi t + \frac{3 \pi}{7})$. And considering your correct decomposition of it into complex exponentials:

$$x(t) = \frac{5 e^{i \frac{3 \pi}{7}}}{2} e^{i 800\pi t} + \frac{5 e^{-i \frac{3 \pi}{7}}}{2} e^{-i 800\pi t}$$

You next apply the continuous-Fourier transform properties to get:

$$X(\Omega) =  5\pi \left( e^{i \frac{3 \pi}{7}} \delta(\Omega - 800\pi) + e^{-i \frac{3 \pi}{7}} \delta(\Omega + 800\pi) \right)$$

Now your Fourier transform has two impulses with weights of $5\pi e^{i \frac{3 \pi}{7}}$ and $5\pi e^{-i \frac{3 \pi}{7}}$ respectively.

The spectral amplitudes are infinite, by definition of the impulse function, but they have the weight respectively. And for the magnitude you would take the absolute value of the amplitudes, which is therefore infinite. Yet again you can represent the magnitudes of those impulses by their weights' absolute value.


How does Music algorithm differ from simple DFT/FFT for frequency estimation?

So both of MUSIC and DFTs are different tools that can be used to accomplish similar goals. DFTs lean on the solid background of Fourier analysis, whereas iterative methods like MUSIC and similar sub-space/AR methods try for better resolution than a DFT is capable of.
With a DFT, the spacing of the frequency bins is prescribed by the number of samples you have, which means that you’ll have spaces “between” your frequency bins that will be ambiguous/smeared between two bins. With MUSIC and similar algorithms, we search over an arbitrary grid of frequencies and can get as fine of resolution as we’d like. In this sense, people will often say that MUSIC is a “super resolution technique”, as it can better resolve closely spaced sinusoids that a DFT would not reveal as being two separate signals.
MUSIC generally outperforms a DFT when trying to detect sinusoids and you know the number of signals a priori. Of course, the major downside here is that you need to tell the algorithm how many signals you care about are present in your measurement. If you don’t know this a priori, it can be difficult to get optimal results.
You should give the Wikipedia page for the MUSIC algorithm a read if you haven’t already, it does a great job comparing the method to the DFT.


Fourier Synthesis With Varying Frequencies Over Time

You might want to have a look at the overlap-add (http://en.wikipedia.org/wiki/Overlap%E2%80%93add_method) and overlap-save (http://en.wikipedia.org/wiki/Overlap%E2%80%93save_method) methods. But, if all you are trying to do is additive synthesis, you don't absolutely have to use the IFFT to generate your signal.

You can set up a bank of $M \in {[1 \ldots N_{Osc}]}$ elementary oscilators of the form $x_m(n)=\alpha_m \times sin(\frac{2 \times \pi \times f_m \times n}{Fs}+\phi_m) $ where $m \in M, N_{Osc} \in \mathbb{N}$, $\alpha_m, \phi_m, F_s, f_m \in \mathbb{R}$ and final output $y(n) = \sum_{m \in M}{x_m(n)}$.

In this case, $\alpha_m, \phi_m$ are amplitude and phase coefficients for each oscillator at frequency $f_m$ and to make them variable in time, all you have to do is turn them to $\alpha_m(n), \phi_m(n)$, i.e. to make them dependent on time sample $n$. This can be an absolute definition, by actually defining a varying waveform or by using some form of interpolation to provide initial and final values for these quantities and let the intermediate values be calculated automatically.

This would work both for generating one sample at a time or a frame of samples at a time (by iteration) and would not have any problems with continuity.


When sinusoidal input starts at n=0, why are transient response associated with z-transform poles of digital filter?

Yes, that's true if the $\mathcal{Z}$-transform of the input signal has poles on the unit circle (i.e., the input is a constant or a sinusoid, starting at $n=0$), and if the filter is stable, i.e. all its poles are inside the unit circle. If you split up $Y(z)$ by a partial fraction expansion, you get a contribution for each pole. The contributions from the poles of $H(z)$ all decay to zero because we've assumed that the filter is stable, i.e. all its pole magnitudes are smaller than $1$, and consequently, these terms decay exponentially. On the other hand, the contributions from the poles of $X(z)$ do not decay because the corresponding poles lie on the unit circle. This means that after a while, all contributions associated with the poles of $H(z)$ will be very small (practically zero; transients), whereas the contributions from the poles of $X(z)$ remain at a constant amplitude (steady-state response).


Reduce "flicker" in spectrum analysis / spectrogram

My suggestion is to low pass filter the FFT slices.

Suppose your current FFTs are:
$$
X(n,k) = {\tt FFT}\{x[n] \cdot w_n[n]\}
$$
where your signal is $x[n]$ and the window $w_n[n]$ enforces how much of the signal you want to include in the FFT at time $n$.

Then form:
$$
Y(n,k) = \alpha \cdot Y(n-1,k)  + (1-\alpha) \cdot X(n,k)
$$
where $\alpha$ is a fraction between 0 and 1 (probably around 0.7 to 0.99).

This will effectively apply a first order lowpass filter to your FFTs.

Another possibility:

$$
Y(n,k) = {\tt median}\{X(n,k)\}
$$

where the median is taken over the last $M$ values of $X(n,k)$ i.e. $\{ X(n,k), X(n-1,k), \ldots, X(N-M+1, k)\}$.


Is it a correct way to calculate the Mean Square Displacement as function of time?

According to the article on which you base your MSD computation, it's specified as a single number (an average over N items) which is defined as the mean (average) of the sum of squared-distances between positions of $N$ particles at time $t$ and given reference positions $x_n(0)$ for each particle at a reference time.

Assuming that you want to compute MSD over a 1D particle displacement of $N$ items, you would use the following Matlab/Octave code to achieve it:

N = 1024;                 % Number of particles
std1 = 2.34;              % set some variance
xpos = std1 *sqrt(3)*2*( randn(1,N) - 0.5); % N particles with positions.
xref = 0;                 % using same reference position for each particle.
MSD = (1/N)* sum ( (xpos-xref).^2 );  % compute MSD of N particles at time t

figure,plot(xpos);        % plot positions of each particle at time t
hold on;                  % plot MSD of those N particles at time t.
plot(MSD*ones(1,N));      % MSD*ones(1,N) is used to plot N points for MSD
title('MSD over particle positions at time t');

Assuming that you also want to plot MSD as a function of sampled time, you could then then use the following code to achieve it:

N = 1024;                 % Number of particles
M = 128;                  % Number of MSD computations sample time points
std1 = 2.34;              % Set some variance, this can change over time.
xpos = std1 *sqrt(3)*2*( randn(N,M) - 0.5); % N particles for each M time samples.
xref = 0;                 % using same reference position for all time and each poarticle.
MSD = (1/N)* sum ( (xpos-xref).^2 );  % now MSD is an 1 x M vector

figure,plot(MSD);        % plot particle positions and MSD


Can't get FFT of Sample & Hold sinusoidal wave in Simulink

I'm no Simulink expert, but I'm pretty sure you don't want sample and hold, but just sample as functionality: As far as I can guess, sample and hold doesn't convert from simulatedly continuous signal to simulatedly discrete signal, but just actually does what its name says, and probably generates a continuous signal, not a discrete one.

In fact, seeing that you're already using a DSP sine, I'd say: instead of using a separate block with a trigger, just change that DSP Sine's sampling rate.


How to Calculate Local Spatial Variance of an Image?

Usually it means the variance within a window of pixels where the (i, j) pixel is in its center:

In the above you can see a 3 x 3 window (Radius of 1).
You'd calculate the variance in this window and set its value to the location of the Pij pixel.


I am currently implementing acoustic FSK modulation and demodulation. I am not a signal processing guy…

Since you say you have matched filters, and you mention non-coherent detection, I think you're pretty much of a digital communication person already – the step to being a DSP person is pretty small :)
The fully-fledged synchronizer SDR approach
So, the classical way to do this is to take the outputs of your two filters, (magnitude) square them, low pass filter the result to something above $2 f_\text{symbol,max}$, and during/after filtering decimate (if possible), so that we don't waste operations per second.
I made myself a quick real-value bandpass signal FSK demod (I think if I did this again, I would have replaced the two bandpass filters with one combined frequency translation and one low pass filter, giving me complex baseband, or threw a polyphase filterbank at the problem, but whatever) in GNU Radio with the help of the GNU Radio companion (flow graph file here):

Point is that we probably want to symbol timing synchronize between the divide and the threshold block. We could do afterwards – the classical, "microcontroller & digital hardware design" approach would in fact be Muller & Mueller clock sink and that would work there just as well – but let's not waste the slope information.
So, I'd add a polyphase clock sync here; I must admit I could not formulate this better than Tom, a disciple of harris, did in the documentation of the eponymous GNU Radio block:

Timing synchronizer using polyphase filterbanks.
This block performs timing synchronization for PAM signals by minimizing the derivative of the filtered signal, which in turn maximizes the SNR and minimizes ISI.
This approach works by setting up two filterbanks; one filterbank contains the signal's pulse shaping matched filter (such as a root raised cosine filter), where each branch of the filterbank contains a different phase of the filter. The second filterbank contains the derivatives of the filters in the first filterbank. Thinking of this in the time domain, the first filterbank contains filters that have a sinc shape to them. We want to align the output signal to be sampled at exactly the peak of the sinc shape. The derivative of the sinc contains a zero at the maximum point of the $\text{sinc}$ ($\text{sinc}(0) = 1,\, \text{sinc}'(0) = 0$). Furthermore, the region around the zero point is relatively linear. We make use of this fact to generate the error signal.
If the signal out of the derivative filters is $d_i[n]$ for the $i$th filter, and the output of the matched filter is $x_i[n]$, we calculate the error as:
$$e[n] = \frac{\Re\{x_i[n]\} \cdot \Re\{d_i[n]\} + \Im\{x_i[n]\} \cdot \Im\{d_i[n]\}}{2.0}\text{ .}$$
This equation averages the error in the real and imaginary parts. There are two reasons we multiply by the signal itself. First, if the symbol could be positive or negative going, but we want the error term to always tell us to go in the same direction depending on which side of the zero point we are on. The sign of $x_i[n]$ adjusts the error term to do this. Second, the magnitude of $x_i[n]$ scales the error term depending on the symbol's amplitude, so larger signals give us a stronger error term because we have more confidence in that symbol's value. Using the magnitude of $x_i[n]$ instead of just the sign is especially good for signals with low SNR.
The error signal, $e[n]$, gives us a value proportional to how far away from the zero point we are in the derivative signal. We want to drive this value to zero, so we set up a second order loop. We have two variables for this loop; $d_k$ is the filter number in the filterbank we are on and $d_\text{rate}$ is the rate which we travel through the filters in the steady state. That is, due to the natural clock differences between the transmitter and receiver, $d_\text{rate}$ represents that difference and would traverse the filter phase paths to keep the receiver locked. Thinking of this as a second-order PLL, the $d_\text{rate}$ is the frequency and $d_k$ is the phase. So we update $d_\text{rate}$ and $d_k$ using the standard loop equations based on two error signals, $d_\alpha$ and $d_\beta$. We have these two values set based on each other for a critically damped system, so in the block constructor, we just ask for "gain," which is $d_\alpha$ while $d_\beta$ is equal to $\frac{\text{gain}^2}{4}$.

Realizing this is audio frequency
So, parameterizing this block properly (and probably reducing bpf_decim andlpf_decim to give the synchronizer more headroom to shift phases) you could build a very stable symbol timing recovery, that would be totally overkill for your application
Since you're doing sampling with a sound card on both ends, and since I presume you don't have to deal with Doppler, the symbol rate has a fixed relation to the carrier frequencies (e.g. the transmitting side has a higher carrier with a period of 8 sound card samples, and a symbol always takes let's say 128 samples, so there's a fixed ratio), you can take a short cut:
If you know the ratio between carrier frequency and symbol rate, just estimate one and you get the other for free! In this case, FSK, that's particularly easy: just wait for a period when the e.g. upper band is active, and use a quadrature demodulator (i.e. practically $\frac{d\arctan}{dt}[n]$) to give you an estimate of the carrier frequency after your decimating band pass filter, and use the resulting value as frequency correction factor – you immediately solve the carrier frequency recovery and symbol rate recovery problem at once. I don't even think you'd need timing recovery after that – after a thresholding, if you oversample enough (which, at these rates, isn't a problem), you can just make a majority decision in a window, or shift that window as long as necessary to get a maximum in the window.
Regarding Packet Sync
Well, since you've got a working demodulator that will probably run too much out of sync when there's no signal present, just use a fixed known data sequence to find your packet.


Matlab: plotting magnitude and phase response

@Hal: Unlike MBaz, I realize you weren't asking for information on the theory of signal processing. But, rather, you were asking a Matlab implementation question. I only have an old version of Matlab so I couldn't try Anthony's code. (Anthony was smart to insert that zero-valued coefficient.) But to plot the frequency mag. and phase responses of a discrete version of your filter, try these:

figure(1), freqz([2,0,3.125],[1,-0.9,0.81]);

or

freqz([2,0,3.125],[1,-0.9,0.81], 512, 'whole');


I don't know about opencv's implementation, but a median filter can indeed be very slow because for every windows you have to sort your data in order to find the median. However there is a recursive implementation of the median filter, which is much more quick. It is very simple to build.

The basic idea is to recursively build the data's histogram that will allow you to very quickly get the median without sorting :
For the first window, fully build the histogram, then, shift your windows by one step (or more, it depends on you), remove the last value from the previous window and add the new one, you updated your histogram with only 2 operations.
Then you have to find the median within the histogram, this is recursively feasible but I don't remember the exact implementation (it is a bit more tricky), however simply checking the histogram's cumulative sum, while brutal, will give you the median quickly (worst case scenario is N summation, with N your window's size).


Correct method for processing FFT in realtime

There are numerous ways to do this and you have the general flow except, as noted, complex averaging is only done in special circumstance, which you haven’t mentioned and are in most cases not applicable.

One can average magnitude or magnitude squared. There is a 1975 Bell System Technical Journal article that goes into detail but a simple comparison is good enough for most people.

Zero padding a FFT to a power of 2 is probably unnecessary if you use a modern FFT library like fftw.  A library will also offer an optimized real FFT. The other advantage of not using a zero pad is that you can pick your bin widths more intelligently.  It may be convenient for a user to have something like 1Hz or 5Hz or 10Hz bin widths.

The number of averages should have some reasoning as well.

Depending on your actual application, some filtering prior to your FFTs might be useful. A seismologist might be interested in frequencies below 10Hz, but a lot of other people don’t. Some degree of pre whitening is often useful.

One can have a “correct” processing chain , but nevertheless underachieve as an analysis tool.

There is a lot about how to display this information that hasn’t been touched. Real time usually means real time display.

When I started doing DSP there were a lot of tricks that you can use to save cycles. You can avoid square roots by taking the max( real, imag) - \alpha min(real, imag) where you set any negatives to positives first.

Another trick is to complex bandshift filter to use a complex fft without redundancy


Preamble Detection and Frequency Offset Estimation in radio packets

For maximum sensitivity preamble detection with time referencing, consider using a barker code, or for more choices with longer lengths a PRN (pseudo-random noise) sequence, or for even more choices when more distinction is needed (CDMA) Gold codes are a possibility.  (Also not further detailed but have other advantages are Kasami codes and Zadoff-Chu codes).

All these codes offer the characteristic of having a strong autocorrelation at time offset = 0 and importantly a much smaller autocorrelation for any other time offsets. Such as that demonstrated in the graphic below for a PRN sequence using an LFSR (Linear Feedback Shift Register):

The codes provide some noise immunity through processing gain; under condition of white noise, the processing gain which is the improvement of SNR is given as $10Log_{10}(N)$ where N is the length of the sequence. This comes out from addition of correlated signal components results in the standard deviation of the signal component growing at rate N while the noise (if white, meaning uncorrelated from sample to sample) will have a standard deviation that grows at $\sqrt{N}$.  We detect the presence of the sequences using correlation, which is a multiply and sum:

$$Corr(\tau) = \Sigma r(t-\tau)x(t)$$

With the formula above representing a cross-correlation function of the received signal r(t) with the seequence of interest x(t), at various time offsets $\tau$.

There are several ways this is done in practice and depends on the detailed implementation approach but for shorter sequences this can be done with an FIR filter, with the coefficients of the filter as the code sequence upsampled to match the sample rate of the received signal (best to do this at the lowest rate possible of course).  For longer sequences it is often more practical to create a code generator with a programmable offset (to shift the start position in time) and simply multiply the received sequence by the code and accumulate the output. A threshold detection can indicate the successful presence of the code.

Below shows the result of a sliding correlation on a GPS signal, where the code repeats every 1 ms. (I would never use an FIR filter to correlate a GPS sequence which is 1023 chips long, but this shows what a sliding correlation would look like).

Other approaches involve doing the correlation in the frequency domain using this relationship:

$$Corr_{xy} = IFFT(FFT(x)FFT^*(y))$$

Which states that taking the inverse Fourier Transform of the complex conjugate product of the Fourier Transforms of x and y results in the circular correlation of x an y.

Further, a very important take-away is the sensitivity of any of these correlations to frequency offset, specifically that the magnitude goes down as a Sinc function with the first null at 1/T.

Please see this post directly where I detail that: GPS CA Signal Acquisition

This informs the frequency capture range based on the length of code used. If you are well inside this main lobe in frequency (I typically use 1/2 the null-null frequency range, but if the SNR is poor you may need to be more constrained). Typical acquisition approaches will either slowly ramp the frequency offset or step the frequency offset during correlation (The frequency is offset by using an IQ rotator on the received signal prior to correlation). Stepping the offset involves setting the frequency at an initial guess, correlating and if the acquisition threshold is not exceeded, stepping the frequency by 1/T (if using 1/2 the null-null frequency range) and repeating. If ramping, the ramp rate must be much less than 1/T without inducing an additional frequency offset loss.

The frequency offset is accurately estimated by comparing two successive correlations, since frequency is the change in phase versus the change in time ($d\phi/dt)$ by simply comparing the phase and time between two successive correlations, we can determine the frequency offset assuming we are in the frequency range of acquisition as detailed above (importantly these correlations are complex correlations with two real outputs I and Q such that carrier phase is determined).

Note that the frequency computation detailed at this link Obtaining phase and magnitude of 2 spectral components after ADC
using the computation

$$K = I[n]Q[n-1] - I[n-1]Q[n]$$  where K is propotional to frequency offset equally applies to the successive I and Q correlation outputs, referred to as I[n], Q[n] and I[n-1]Q[n-1].

Much more practical with modern processing and shorter preamble sequences would be a one shot frequency offset and delay estimation that I have detailed here:

GPS signal acquisition

That said, a robust preamble solution could involve sending multiple copies of one code to aid initial carrier and clock acquisition (and channel estimation!) followed by a final negative of the code sequence to mark the beginning of the data. (Just invert the code and the correlation will be shifted 180°). The trade space is over-head of the preamble versus acquisition sensitivity. It is typical to require more SNR during acquisition, especially if a fast acquisition is needed, hence the significance of correlation and processing gain.


Compressive Sensing - Location of Non Zero Elements

Actually some greedy algorithms like MP, exactly do that. They iteratively find signal supports (location of non-zero elements) using measurements and measurement matrix, and reconstruct original signal.
However, in case you are asking for a method to find signal supports from measurements all at once, I think there is no such thing (may be someone discover it someday!).
If there was a such a thing, the computational cost of the CS recovery algorithms wouldn't be so high.

Fore further information on greedy recovery algorithms, take a look at these:

http://www.mathworks.com/help/wavelet/ug/matching-pursuit-algorithms.html

https://en.wikipedia.org/wiki/Matching_pursuit


removal of blackman harris window results in signal higher than expected

The reason for that is that you don't normalize the DFT samples properly. Dividing by number of samples in time domain is valid only for Rectangular Window. For simple case of DFT, you should divide your amplitude by sum of window samples:

$$S=\sum_{i=0}^{N-1}w_i $$

This for Rectangular Window will be obviously equal to $N$ ($8192$ in your case), whereas for Blackman-Harris Window it's about $2939$, therefore the Coherent Gain is about $0.36$ instead of $1$.

For reading more about this interesting topic I recommend you these papers (especially the first one):

  G. Heinzel et al - Spectrum and spectral density estimation by the
  DFT, including a comprehensive list of window functions and some new
  flat-top windows - especially chapter 9

  F. J. Harris - On the Use of Windows for Harmonic Analysis with the
  Discrete Fourier
  Transform


Discrete-time Fourier Transform of the unit step sequence $u[n]$

Cedron Dawg posted an interesting initial point in this answer. It begins with these steps:

$$ \begin{align}
U(\omega) &= \sum\limits_{n=0}^{+\infty}  e^{-j \omega n} \\
&= \lim_{ N \to \infty } \sum\limits_{n=0}^{N-1}  e^{-j \omega n}\\
&= \lim_{ N \to \infty } \left[  \frac{  1 -  e^{-j \omega N}  }{ 1 -  e^{-j \omega }   } \right] \\
&= \frac{ 1 }{ 1 -  e^{-j \omega } } - \lim_{ N \to \infty } \left[  \frac{  e^{-j \omega N}  }{ 1 -  e^{-j \omega } } \right]
\end{align} $$

It turns out the term inside the limit can be expanded as follows:

$$\begin{aligned}
\frac{  e^{-j \omega N}  }{ 1 -  e^{-j \omega } }
={}
&\frac{1}{\sin^2(\omega)+(1-\cos(\omega))^2} \cdot \\
&\left[-\cos(\omega)\cos(N\omega)+\cos(N\omega)-\sin(\omega)\sin(N\omega)+ \\
j(-\sin(\omega)\cos(N\omega)+\cos(\omega)\sin(\omega)-\sin(N\omega))\right]
\end{aligned}
$$

The common factor outside the brackets can be expressed as:

$$\frac{1}{\sin^2(\omega)+(1-\cos(\omega))^2} = \frac{1}{4\sin^2(\omega/2)}$$

The real part inside the brackets also equals:

$$-\cos(\omega)\cos(N\omega)+\cos(N\omega)-\sin(\omega)\sin(N\omega)= 2\sin(\omega/2)\sin[\omega(-N+1/2)]$$

On the other hand, the imaginary part can be rewritten as:

$$-\sin(\omega)\cos(N\omega)+\cos(\omega)\sin(\omega)-\sin(N\omega)=-2\sin(\omega/2)\cos[\omega(-N+1/2)]$$

Rewritting the original term we get that:

$$\begin{align}
\frac{  e^{-j \omega N}  }{ 1 -  e^{-j \omega } }
&=\frac{2\sin\left(\frac \omega 2\right)}{4\sin^2\left(\frac \omega 2\right) }
\left( \sin[\omega(-N+1/2)] - j\cos[\omega(-N+1/2)]\right) \\
&=-\frac{\sin[\omega(M+1/2)]}{2\sin\left(\frac \omega 2\right)} - j\frac{\cos[\omega(M+1/2)]}{2\sin\left(\frac \omega 2\right)}
\end{align}
$$

where I used $M=N-1$ and the limit stays unaffected as $M\to \infty$ as well.

According to the 7th definition in this site:

$$\lim_{M\to \infty} -\frac{1}{2\sin(\omega/2)}\sin[\omega(M+1/2)] = -\pi \delta(\omega)$$

So far we have that:

$$\lim_{ M \to \infty }   \frac{  e^{-j \omega (M+1)}  }{ 1 -  e^{-j \omega } } =-\pi\delta(\omega)
-j\lim_{M\to \infty} \frac{\cos[\omega(M+1/2)]}{2\sin(\omega/2)}
$$

If we could prove that the second term on the right of the equality is $0$ in some sense, then we are done. I asked it at math.SE and, indeed, that sequence of functions tends to the zero distribution. So, we have that:

$$ \begin{align}
U(\omega) &= \frac{ 1 }{ 1 -  e^{-j \omega } } - \lim_{ N \to \infty } \left[  \frac{  e^{-j \omega N}  }{ 1 -  e^{-j \omega } } \right] \\
&=\frac{ 1 }{ 1 -  e^{-j \omega } }+\pi\delta(\omega)
+j\lim_{M\to \infty} \frac{\cos[\omega(M+1/2)]}{2\sin(\omega/2)} \\
& =\frac{ 1 }{ 1 -  e^{-j \omega } }+\pi\delta(\omega)
\end{align} $$


Fourier Transform of Kernel Density Estimation - Convolution Theorem?

It's a bit contrived, but observe that by the "sifting property" of Dirac-delta function:

$$
K(x-X_j) = K(x) * \delta(x-X_j)
$$

where $*$ denotes convolution and $\delta$ is the Dirac-delta function.

Now you can apply Fourier convolution theorem and write:

$$
\mathcal{F}\{K(x-X_j)\}(t) = \kappa(t) e^{j t X_j}
$$

where $\mathcal{F}\{g(x)\}(t)$ denotes the Fourier transform of a function $g(x)$.

Finally using the linearity property of the Fourier transform:

$$
\mathcal{F}\{\hat{f}(x)\}(t) = \frac{1}{N} \sum_{j=1}^N \mathcal{F}\{K(x-X_j)\}(t) = \kappa(t) \sum_{j=1}^N e^{j t X_j}.
$$


Amplitude demodulation of sinusidal signal in matlab

I think there are some problems with your analysis structure.

First of all, for being able to see from -0.5 to 0.5, without aliasing, you have to sample at 1.

Secondly, if you want to down convert a signal, bring it to DC you need to multiply it with the corresponding cosine in time, or discrete time, or apply the convolution to its spectrum with the corresponding Fourier transformed cosine (a cosine in spectrum is a delta Dirac function on the frequency of the same cosine).

By the way, I do not know what the mkhwdata function does.
The solution should look something like this:

clc
close all;
clear

fc = 0.2;
fs = 1; %%%%%%%%%%%%%%%%%%

[z]= mkhwdata([xxxxxxxxx]);

[Z,nu] = tdftfast(z);
%N = length(z);
n = [0:(length(z) - 1)];
%n = N*(0:N-1);
mult  = exp(1i*2*0.25*pi.*n);
mult2 = exp(1i*2*pi*0.5.*n);

% We bring the image to DC
[z1] = (z.*transpose(mult));
[z2] = (Z.*transpose(mult2));

%Lowpass filter
[num,den] = butter(9,0.2);

y  = filter(num,den,z);
y1 = filter(num,den,z1);
y2 = filter(num,den,z2);

figure(4)
ylabel('Amplitude');
xlabel('Frequency (hz)');
subplot(221)
present_image(z);
subplot(222)
present_image(y);
subplot(223)
present_image(y1);
subplot(224)
present_image(y2);


result = [math.degrees(x/y)*(-1) for x,y in zip(BxHPF,ByHPF)]

It seems that I have found an answer. Maybe it will help someone.


Confusion about ensembles and averages in autocorrelation matrices

EDIT: sorry that I didn't fully get your question. Now I got your points of conufion and try to provide a shortest possible answer:

First, you should better use xcorr function when estimating ACS $r_x[k]$ from data sequence $x[n]$. The output of xcorr will return the lag-0 at the middle sample as you recognize. Look at the input argument scaleopt which provides normalization to the computed output. The default is none ! and causes a confusion for you. You should select biased so that equation 7 below and xcorr output will be identical. Like this:

rx = xcorr([1,2,3],'biased'); % identical to (1/3)*xcorr([1,2,3])

Second,  do not use the raw output of xcorr as the input to the toeplitz function. As the output of xcorr will include lags for rx[-2],rx[-1],rx[0],rx[1],rx[2]. But you only need rx[0],rx[1],rx[2]. So you better call it like:

L = 3;
x = [1,2,3];
rxL = xcorr( x, 'biased');
Rx = toeplitz( rxL(L:end) );  % now yo uare ok ;-)

Now it could be argued that the Matlab's default xcorr scale option none is misleading. And yes, I would say so. In its default form it computes the convolution and leaves it unscaled...

Finally, please distingusih a theoretical expectation from practical estimation. When algorithms work on data, they estimate and they don't compute expectations, which is a theoretical tool. Hence, all expectations are replaced with estimations when algorithms are run for the practical applications...
It seems (to me) that your confusion is about three things.

Theoretical definition of the autocorrelation sequence $r_{XX}[k]$ of a random process $X[n,s]$
Practical estimation of it from a sample function $x[n]$ via ergodicity.
Interpretation of a random vector (including samples of a random process)

Consider the interpretation of a random process $X[n,s]$ as an ordered collection (set) of indexed random variables $X_n(s)$ for each time index $n$. Each random variable $X_n$ may, in principle, has its own unique probability distribution/density function $F_{X_n}(x)$. Consider the enormous implication of this on the number and type of all possible joint density functions.

To make things extremely simplified, we make the assumptions that all those random variables have the same PDF $F_{X_n}(x)$ for all $n$, and independent of each other, also further that first and second moments are independent of partciular time instants, which yields an i.i.d. , WSS random process $X[n,s]$ at hand.

The theoretical auto-correlation sequence (ACS) of this WSS random process is defined as:

$$ r_{X_n X_m}(n,m) = E\{X_n X^*_m\} = E\{X[n]X^*[n+k] \} = r_{X X}[k]   \tag{1}$$

The computation of this theoretical ACS $r_x[k]$ requires you to evaluate the integrals with PDF functions. This is possible with questions on the paper, but with applications in working systems, where you won't solve (mathematical) questions but perform statistical calculations, you need to estimate it from data, based on ergodicty assumption, which states that theoretical ensemble averages can bu estimates from practical sample function time-averages. Then you have

$$ \hat{r}_X[k] = \lim_{N \to \infty} \frac{1}{2N+1} \sum_{n=-N}^{N} x[n]x^*[n+k] \tag{2}$$

Note that $x[n]$ is a particular sample function (a realization) of the hrandom process $X[n,s]$. Note also that the infinite limits are replaced by data observation boundaries in practice.

Now consider the random vector $X(n) = [X(n),X(n-1),..,X(n-N+1)]^T$ which includes $N$ consecutive random variables from the random process $X[n,s]$, beginning at index $n$.

From this vector build up the outer product $X \cdot X^H$ , which includes a matrix of pairwise multiplications of random variables $X_i X_j$ as in

$$
RM =
  \begin{bmatrix}
X_1 X^*_1 & X_1 X^*_2 & ... & X_1 X^*_N \\
X_2 X^*_1 & X_2 X^*_2 & ... & X_2 X^*_N \\
... & ... & ... & ... \\
X_N X^*_1 & X_N X^*_2 & ... & X_N X^*_N \\
  \end{bmatrix} \tag{3}
$$

Now taking the expectation of this $RM$ matrix yields the theoretical auto-correlation matrix of the WSS random process as:

$$ R_X = E\{ X \cdot X^H \} = E\{RM\} = E\{\begin{bmatrix}
X_1 X^*_1 & X_1 X^*_2 & ... & X_1 X^*_N \\
X_2 X^*_1 & X_2 X^*_2 & ... & X_2 X^*_N \\
... & ... & ... & ... \\
X_N X^*_1 & X_N X^*_2 & ... & X_N X^*_N \\
  \end{bmatrix} \} \tag{4}$$

and this results in:
$$ R_X = E\{ X \cdot X^H \} = \begin{bmatrix}
r_x[0] & r_x[-1] & ... & r_x[-N+1] \\
r_x[1] & r_x[0] & ... & r_x[-N+2] \\
... & ... & ... & ... \\
r_x[N-1] & r_x[N-2] & ... & r_x[0] \\
  \end{bmatrix} \tag{5}
 $$

for real data it becomes:
$$ R _x =
 \begin{bmatrix}
r_x[0] & r_x[1]^* & ... & r_x[N-1]^* \\
r_x[1] & r_x[0] & ... & r_x[N-2]^* \\
... & ... & ... & ... \\
r_x[N-1] & r_x[N-2] & ... & r_x[0] \\
  \end{bmatrix} \tag{6}
$$

Again this theoretical ACM $R_X$ cannot be computed but is estimated from sample function data...

Estimation of the auto-correlation matrix $R_x$ can be done in a number of ways (both related) but at the basic level you see that what you should actually estimate is the auto-correlations sequence $r_x[k]$ for the lags $k=0,1,...,N-1$. Then from the Toeplitz and Hermitian symmetric property of $R_x$, you can construct it.

One way to practically estimate ACS $r_x[k]$ from the data $x[n]$ of length $L$ is :

$$ \hat{r}[k] = \frac{1}{L} \left( x[k] \star x^*[-k] \right) \tag{7}$$

This yields an estimate of ACS for the lags $-L < k  < L$, (you only need $N$ of those for $k=0,1,...,N-1$)

Alternatively, ACM can be estimated as a whole similarly to ACS from:

$$ \hat{R}_x = \frac{1}{K} \sum_{k=1}^{K} X(k) \cdot X(k)^H \tag{8}$$  where data vector $X(k)$ includes $N$ consecutive samples from $x[n]$ (specially processed at the edges)


How to narrow the frequency band of a wavelet

Center frequency and bandwidth are linked in a single "scale" parameter for wavelets. That means if you change the bandwidth you also change the center frequency and vice versa.

For the discrete wavelet transform on a dyadic lattice the relative bandwidth is also fixed at approximately one octave.

If you want a continuous wavelet transform frame with a smaller relative bandwidth you need to construct a different wavelet. Intuitively, it will come with more oscillatory cycles, i.e. more zero crossings.

The Ricker wavelet you have coincides with the 2nd so called Hermite function. Higher order Hermite functions also have an increasing number of zero crossings. Unfortunately they're one of the examples of where the intuition of adding more cycles will not work. The Hermite functions grow in both frequency as well as time support with increasing order. So that's not the way to go. I'm only making this clear because you might find a reference to them in relation to the Ricker Wavelet.

So instead I would try to adapt the construction of the wavelet to your specific needs, which you will have to outline in your question to receive a more detailed answer.


How to make sure that a sinusoid wave with different subsequent frequencies matches up?

You need to make sure that the phase of continuous.

The easiest way would be something like this.

x(n) = cos(phase);
phase = phase + 2*pi*current_frequency*sample_time;
current_frequency= update_frequency(n);
n = n + 1;

Note, that this is very inefficient code and only for illustrative purposes.


Scipy Butterworth highpass filter results in distortion and bloated output

For digital filters, Wn are in the same units as fs. By default, fs is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (Wn is thus in half-cycles / sample.)
For analog filters, Wn is an angular frequency (e.g. rad/s).
As yours is an analog filter, your Wn parameter is 2π·54Hz
from scipy import constants
...
Wn = 2*constants.pi*54
b, a = signal.butter(order, Wn, btype='highpass', analog=True)


Short answer is no, there are none, or maybe none that made it until becoming a de facto standard like FAST or SIFT.
This is still an open area of research by the way if you plan to work in this field.

Mostly, the issues are linked to the fragility of color representation in digital images. The common RGB color space is not really relevant (distances between colors in RGB do not match the perceptual intuition of color distance), so you probably have to go Lab or YCbCr first, then work on the luminance channel. Also, be aware of the white balance that can vary a lot between different cameras, or with a given camera but varying light conditions (sunlight, artificial light...).

Furthermore, in most cases the detection and/or description have good reasons to be color agnostic. For example, when you watch a black-and-white movie, your brain can still recognize objects although they have lost their color.


conversion of 1d signal to 2d signal for CNN operation

There are a couple things with that question.

1) Why do you need 2D convolutions exactly? You could think of a 1D signal as an $Nx1$ 2d signal if you like ;) Then you can design your convolutional filters to be $7x1$, $5x1$, $3x1$ etc. The pooling you only apply in the x-direction. And there you go. It's a 2D convolutional neural network, which works on 1D signal. Just kidding, it just a CNN.

2) Simple fully-connected deep learning neural networks are better for handling 1D data, as far as I know. Of course it's best if your data is evenly sampled. CNN comes handy when the fully connected structure explodes the computation. In your case, that should not be a big issue.

If your 1d data vector is too large, just try subsampling, instead of a convolutional architecture.

3) Converting a 1d data to 2d is probably valid only if you know in advance that this 1d manifold carries non-uniform neighborhood information, which could be represented with a 2D matrix with nearby connections. The paper probably inserts a layer to the network which re-arranges the indices into a 32x32 matrix. The weights or the indices where to put each value maybe learnt by the network itself. You would need to write the backprop for this though.


How to know when a sound in a given frequency starts and ends?

I would recommend reading about onset detection. It's about music signals but nevertheless a few ideas might be extracted from it.

From a more general point of view, what I would do is a frame-based analysis of the signal and define some kind of energy threshold/noise floor with which I could determine if there is activity or not.


How to Filter Human Speech from an Audio clip?

This is not a complete answer, but detection of words from the speech itself is almost always dictionary based - hence one needs to know what words would sound like. If you don't know language, you don't quite know the relation to the sounds. Hence, speech/audio recognition don't quite work.

There can be many primary features to classify:
Reference:
Automatic Language Identification: A Review/Tutorial  by Yeshwant K. Muthusamy , Etienne Barnard , Ronald A. Cole Signal Processing Magazine, IEEE Oct. 1994 Volume: 11,   Issue: 4 pp. 33 - 41

Or related paper of same author.

Phonetics:
Here you should try to look for is that you should try to collect and identify sounds based on primary phonetics. For example - you can use something like International Phonetic Alphabet
Prosodic:
Languages like English, French, German might have very close phonetics. In such case, you can also use Prosodic information which is essentially duration of sound, emphasis or accentuations, pitch and pitch variations.
Phonotactics:
Phonotactics refers to the rules that govern the combinations of the di erent
phones in a language. There is a wide variance in phonotactic rules across languages. For example, the phone cluster /sr/ is very common in the Dravidian language Tamil, whereas it is not a legal cluster in English.
Vocabulary:
If you are successful at this, you can probably make another classification of which phonetic alphabet orders are most likely in particular language. You probably don't have to recognize each word. But based on most likely phonetic symbol sequences, you can probably do language detection.


It is going to be a tough task, and you may not achieve an automated process.

Here are some bricks that you could try:

Enhance the contrast of the image so that you can clearly distinguish all the features by looking at the image. Here is a link to an algorithm that shows good results.
In this enhanced image, extract the edges (using for example a standard Canny-Deriche edge detector).
From the edge detector output, you can start adding the required application logic in order to keep the desired vessels: size or length criteria, linear features extraction, etc.


Would it be feasible to implement the Sliding DFT on a digital signal processor?

The paper Sliding DFT control algorithm for three-phase active power filter implements a sliding a sliding DFT on a TMS320F2812 from Texas Instruments, and discusses some advantages, and An Enhanced Interpolated-Modulated Sliding DFT for High Reporting Rate PMUs does similar work on FPGA.

From what I remember, you may have a look at the stability, see  The Sliding DFT.


You can look at implementations for open source image editing programs (e.g. GIMP). This link from MSDN Magazine, January 2005 by Stephen Toub, entitled "Sepia Tone, StringLogicalComparer, and More" in the .NET matters section specifies a linear transformation of the color space to get a Sepia effect: "New Red = R*.393 + G*.769 + B*.189 New Green = R*.349 + G*.686 + B*.168 New Blue = R*.272 + G*.534 + B*.131" where each channel is normalized to lie in $[0,1]$.

Its just applying the formula on each pixel, so you don't need any special libraries.


Implementation of IIR Filter in dsPIC 33EP - Q15

There's an old trick you can use:

$$1\cdot x=x$$

So you don't need to multiply with a factor that equals $1$.


Applications or physical interpretation of auto-convolution?

Autoconvolution is used in signal detection, but the way you've written it is not correct. Suppose you're trying to detect a signal $f(t)$ by filtering with h(t).

$y(t) = f(t) \ast h(t)$

You want to maximize your response to the signal $f(t)$. We can do this by maximizing the correlation coefficient between $f$ and $h$. Here the correlation is time-varying , so we'll maximize the average autocorrelation coefficient. We'll assume the $f$ and $h$ signals have DC values of zero for simplicity.I'll use $\mu_y$ to denote average value of a the auto-correlation of response, $y$.

$ h = argmax_{h} \ \   \mu_y\big(E_f E_h \big)^{-\frac{1}{2}}$

Let's take a differential of our performance $ J \propto \mu_y\big(E_f E_h \big)^{-\frac{1}{2}}$ with respect to h

$\partial_{h(\tau)} J \propto \partial_{h(\tau)} \bigg((\mu_y\big(E_f E_h \big)^{-\frac{1}{2}}\bigg)$

$\ \ \ \ \ \ \ \ \  =  \frac{1}{|dom(f)|}\bigg(\partial_{h(\tau)}\mu_y\bigg)\big(E_fE_h)^{-\frac{1}{2}} +\frac{1}{2}\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\bigg(\partial_{h(\tau)}E_h\bigg)$

$\ \ \ \ \ \ \ \ \  \propto  \bigg(\partial_{h(\tau)}\frac{1}{|dom(f)|}\int_{\infty}^{\infty}f(t-\tau)h(\tau)d\tau \bigg)\big(E_fE_h\big)^{-\frac{1}{2}} +\frac{1}{2}\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\bigg(\partial_{h(\tau)}E_h\bigg)$

$\ \ \ \ \ \ \ \ \  =  \bigg(\frac{1}{|dom(f)|}\int_{\infty}^{\infty} f(t-\tau)d\tau \bigg)\big(E_fE_h\big)^{-\frac{1}{2}} -\frac{1}{2}\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\bigg(\partial_{h(\tau)}E_h\bigg)$

$\ \ \ \ \ \ \ \ \  =  \big(E_fE_h\big)^{-\frac{1}{2}} \frac{1}{|dom(f)|}\int_{\infty}^{\infty} f(t-\tau)d\tau-\frac{1}{2}\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\bigg(\partial_{h(\tau)}\int_{\infty}^{\infty} h^2(\tau)d\tau \bigg)$

$\ \ \ \ \ \ \ \ \  =  \frac{1}{|dom(f)|}\big(E_fE_h\big)^{-\frac{1}{2}} \int_{\infty}^{\infty} f(t-\tau)d\tau-\mu_yE_f^{-\frac{1}{2}}E_h^{-\frac{3}{2}}\int_{\infty}^{\infty} h(\tau)d\tau$

$\ \ \ \ \ \ \ \ \  \propto  \frac{1}{|dom(f)|}\int_{\infty}^{\infty} f(t-\tau)d\tau-\mu_yE_h^{-1}\int_{\infty}^{\infty} h(\tau)d\tau$

$\ \ \ \ \ \ \ \ \  = \int_{\infty}^{\infty} \bigg(\frac{1}{|dom(f)|}f(t-\tau)-h(\tau)\mu_yE_h^{-1}\bigg)d\tau$

At minimum $J$ we shouldn't assume $h(\tau) = 0$, we'll enforce $$\frac{1}{|dom(f)|}f(t-\tau)-h(\tau)\mu_yE_h^{-1} = 0$$

and pretty readily you get

$$h(\tau) = \frac{E_h}{\mu_y|dom(f)|}f(t-\tau)$$

Or most importantly

$$h(\tau) \propto f(t-\tau) \text{where} \ \ \tau \ \ \text{is time and} \ \ t \ \ \text{is the delay}$$

or, using perhaps better "variable names"

$$h(t) \propto f(t_d-t) \text{where} \ \ t \ \ \text{is time and} \ \ t_d \ \ \text{is the delay}$$

That is, if we want to maximize the correlation between our signal detector with impulse h(t), we better pick h(t) to be a time-reversed and time-shifted version of our signal of interest. In practice $t_d$ would probably be set to zero, as it just represents whenever the $f$ part that you're looking for finally arrives.

Under this chosen $h(t) \propto f(t_d - t)$, your original question makes more sense. The autocorrelation signal becomes proportional to the accumulated energy of $f(t)$ that is seen by your filter.

$y(t) = f(t) \ast h(t)$

$y(t) \propto \int f(\tau) h(t-\tau)d\tau$

$y(t) = \int f(\tau) f(t_d-(t+\tau))d\tau$

$y(t) = \int f(\tau) f(\tau + t_d - t)d\tau$

$y(t) = \int P_f(\tau + t_d - t)d\tau$


If you can get input signal when it doesn't contain any useful signal (I mean only noise is presented), you can estimate average noise power at first. Simply find a power of such a signal:

$P_n=1/N \cdot \displaystyle\sum_{k=0}^{N-1}|s(k)|^2$.

Choose $N$ in $2^{12}\ldots 2^{15}$ range for example.

Then you can measure a power of $signal+noise$ mixture, $P_s$, by the same way, when input signal is presented at the input. Caution: the measurement of $P_s$ is only valid if the signal exists all the measurement time. SNR could be calculated by:

$SNR=10 \cdot \log_{10} \frac {P_s - P_n}{P_n}$.

If your noise is stationary and white it's possible to measure it once in the setup and use $P_n$ as a constant.

If you can't measure noise power at any reason, you can try to estimate SNR in the frequency domain via FFT. But you have to know your signal's band. The point is the same but you should use FFT frequency indexes instead of sample numbers $k$ while estimating integral power. $P_s$ is now power of $signal+noise$ in the band of interest, while $P_n$ is the noise power in total band ($0...fs$). In this case power spectral density of the noise in your system should be evaluated at first. This method could be far too complex for your task.

Both methods show you pointless result in the case of low SNR in the signal band (not total band in general). Frequency domain method will be preferable if your signal's band is quite less than $f_s$, otherwise time domain method will be good enough. To estimate SNR of the order of 5 dB or lower it's better to know signal and noise power (or power spectral density) a priory. So the task has robust analitical solution.

Hope this helps.


how to handle different durations of audio data?

In this case, usually, Normalization is done.
For example in your training and testing data, you need the same shape, so that you should try something like,
mean = np.mean(X_train_features, axis=0)
std = np.std(X_train_features, axis=0)
X_train_features = (X_train_features - mean)/std

Here X_train_features can be a data frame of spectrogram or mfcc features.
The same thing you can apply for testing features. One thing to remember is the number of columns in shape of training and testing should be the same.
You can also check some Kaggle kernels which are related to Audio processing.


cceps (complex cepstrum) function in Matlab

You should split the signal up into overlapping chunks of 20 to 30 ms. It is certainly valid to take the complex cepstrum on the 4 second signal. But, since the frequency characteristics of speech are constantly varying the complex cepstrum on the 4 second signal will average all of these together and you won't know what happened when. If you take it in each 20 ms window then you can see how the properties of the signal change over time.


You are correct. The following is what you have shown.
Let $x[n]$ be a periodic sequence with period N :
$$ x[n] ~=~ x[n + rN] ~~~~~,~~~\forall r,N \in Z^+ $$
Define $y[n] = x[M n]$, where $M$ is a positive integer; then
$$ y[n] = x[Mn] = x[Mn + rN] = x[M(n + \frac{rN}{M})] = y[n+K] $$
$y[n]$ will be periodic with $K$, if $\frac{r ~N}{M}$ is an integer. Hence, set $r$ to any value that makes $\frac{rN}{M}$ an integer.
If $N,M$ are already coprime, then setting $r=M$ makes $K=N$ as the period of $y[n]$. If they are not already coprime, then make them coprime $M_c,N_c$, and then set $r=M_c$ to make $K = N_c$ as the period of $y[n]$.
Make $N,M$ coprime by dividing them by their greatest common divisor $g = \text{gcd}(N,M)$ :
$$ N_c = N / g ~~~,~~~~ M_c = M / g $$
Then, for the general case, the period of $y[n]$ is:
$$ \boxed{ K = \frac{N}{ \text{gcd}(N,M) } } $$


Will similar enough vectors be "put together" and only treated as 1?

No, this is not how it is done.

  so for what p will P(o|p) be computed?

Your assumption that recognition is performed by considering one by one all possible phone or word sequences and scoring them is wrong. For just a few seconds of speech the number of possible words would be extremely large! This approach - scoring each candidate model among the entire space and picking up the one with the largest score, is actually used only for some small vocabulary / single-word applications (say recognizing several options on a voice menu). I think this confusion stems from the fact that the classic Rabiner tutorial on HMM makes a lot of sense for finite-vocabulary applications, but stays away from all the difficulties/technicalities of continuous applications.

In continuous speech recognition there are several layers on top of the phones models. Phone models are joined together to form words models, and words models are joined together to form a language model. The recognition problem is formulated as searching the path of least cost through this graph.

One could do the same by piecing together HMMs, though... Let's take a simpler example... You have 5 phones and you want to recognize any sequence of them (no words models, no language models - let's consider meaningless speech for this example). You train one left-right HMM per phone, say 3 or 4 states - maybe on clean individual recordings of each phone. To perform continuous recognition, you build a big HMM by adding transition states between the last state of each phone and the first state of each phone. Performing recognition with this new model is done by finding the most likely sequence of states given the sequence of observation vectors - with the Viterbi algorithm. You will know exactly which sequence of states (and thus phones) is the most compatible with the observation. Practically, things are never done this way, and this is why I downplayed in my previous reply the importance of Rabiner's "problem 2": in a true continuous application, the size of the vocabulary is such that considering this giant HMM made of all possibles connections between phones to make words, and words to make sentences would make the naive use of the Viterbi algorithm impossible. We stop reasoning in terms of HMM and we start thinking in terms of FST.


Is there a special case in which wave atoms convert to curvelets?

In their 2007 paper, Wave atoms and sparsity of oscillatory patterns, Demanet and Ying  draw connections between wave atoms and other types of wavelets:

In plain words:

  We introduce “wave atoms” as a variant of 2D wavelet packets obeying
  the parabolic scaling wavelength

with a better spatial localization, meant for texture analysis. Their choice ($\alpha=\beta=1/2$) was a compromise between two types of sparsification (preservation under warping, sparsity of oscillations). Yet, one could built different  collections of such wave packets, with alternative  choices of  $\alpha$ and $\beta$, bridging the gap between the different transform choices.

Yet, their natural redundancy is 2,  with an orthonormal variant, and a complex one with redundancy 4, different from that of the curvelets. So I do no think they can be used as a special case of curvelets.

A tutorial on 2D geometric transformation can be obtain in A Panorama on Multiscale Geometric Representations, Intertwining Spatial, Directional and Frequency Selectivity, 2011.


Why does dot product return zero at non-harmonic frequency

Note that when you compute the DFT, you compute the dot product with a sine term and with a cosine term. In your third example, even though the dot product with the sine term is zero, the dot product with the cosine term isn't, so you do get spectral leakage, as you expected. The dot product with the sine term is zero because you have a whole number of periods plus one half period inside the DFT frame, so the periodized version of the signal is even and, consequently, the imaginary part of its DFT (the dot product with the sine term) vanishes.


Best way to remove unwanted region from image

Looks like you are using Matlab. Try bwareaopen(I, N), where I is the original binary image and N is the estimated size of each unwanted connected region.

You can try edit bwareaopen for more details. Basically the algorithm tries to find the size of connected regions. Connected-component labeling with union-find algorithm is expected to get you there.


Soft question: seeking online course reference for CS

EDIT: There is no shortage of links to go through:
http://www.quora.com/Compressed-Sensing/What-are-some-good-resources-for-learning-about-compressed-sensing

If anyone can suggest a good one from experience I would appreciate it.


How to calculate the channel gain if the channel is a $N$ by $M$ matrix,$N,M>1$

It is important to understand the what the channel matrix $\mathbf{H}$ physically means. The channel matrix is a symmetric matrix that contains the fading coefficients, that is, in row $i$ column $j$ the value ${h}_{ij}$ is the coefficient between the $i^{th}$ transmit antenna and the $j^{th}$ receive antenna. The gain between any transmit/receive antenna pair can be found by $|{h}_{ij}|^2$.


Stable digital filter with leading zeroes

OK, so this is a bit of an intuitive argument so I'm not 100% sure about it, and it's totally sideways to the way I normally think about filters, but:

Let's consider the "feedback kernel" ($a_k ... a_n$ - so, not including $a_0$) as an FIR filter. If the magnitude of that kernel in the frequency domain is $< 1$ including inter-bin peaks, the original filter should be stable.

If so: to construct a filter: one could generate a random set of $a_k ... a_n$, find the maximum freq response of that feedback kernel, and scale such that the peak is less than 1, to get a stable feedback kernel.

To reason about this: let's construct an infinite sequence of sequences, $Y^k$:

$$
Y^0_i = a_0x_i\\
Y^j_i = -\sum_{l=k}^na_lY^{j - 1}_{i - l}
$$

That is, each $Y^j$ is the previous $Y^{j - 1}$ convolved with our feedback kernel (the negation doesn't affect the argument).  I also think it's true that:

$$y_i = \sum_j Y^j_i$$

Let's say $x_i$ is "time-bounded" (zero outside a finite range) - it follows inductively that every $Y^j$ is also time-bounded, because it's a convolution of the previous.

Let's also say that the feedback kernel's transfer function has magnitude $< G$ for all frequencies.  This guarantees that whatever the input, the total energy in each $Y^k$ is strictly less than $G$ times the previous one.  This implies a geometrically-decreasing upper bound on the magnitude of the sample values of $Y^j$ as well.

This geometric progression of bounds means that when we sum up $y_i$ it converges if $G < 1$, and to a value bounded by $K/(1 - G)$ for some $K$.  The time-boundedness means that it decays over time (hand-waving over the formalism here).

Does this hold up?


Interpolating and Decimating CIC filter based on sinc function

The frequency/z domain response of a $N$-Stage CIC filter is, invariably [1]:

$$\begin{align}
H(z) &= \frac{(1-z^{-RM})^N}{(1-z^{-1})^N}\\
&= \left(\sum\limits_{k=0}^{RM-1}{z^{-k}}\right)^N
\end{align}$$

with $R$ being the rate change, and $M$ being the delay length. You can assume $M=1$ (you don't get a special shape if you assume larger $M$, just sharper response).

For $N=1$, you get the "boring" observation that you've built a moving average filter (which, in fact, has sinc-shape in frequency domain). But that's really just due to the fact that that you really did just that – build a recursive MA.

The magnitude response is usually (as in [1]) represented as

$$
\lvert H(f)\rvert = \left\lvert \frac{\sin(\pi M f)}{\pi M f} \right\rvert^N
$$

which goes to say you only get the sinc shape for $N=1$. In other words, a moving average. Which really shouldn't surprise you – the moving average is the time-domain rectangle, and the Fourier transform of that is the sinc.

[1] CIC Filter Introduction, Matthew P. Donadio, July 2000. Available Online


[First shot] What I like with SE.DSP is that you can try to answer questions on topics you know nothing about. And I know nothing about space-time block coding. But I'll dig into the papers to be more precise later, since I am quite interested in concepts like diversity, redundancy, sparsity.

Most signal processing algorithms reside in the gray area between: "what you really want" and "what you can reasonably achieve". Sometimes, what you want is not computable, or difficult to model, and most of your DSP  work consists in finding a "good-enough" proxy to satisfy both goals. This the modelling and optimization part of DSP.

Concepts like diversity, redundancy, sparsity are rather vague, and could formalized in many and very different ways, depending in the objective. One acceptation of diversity is: how good am I at describing the richness of the space of possibilities? Good may means: do I describe the space of possibilities as much as I can (with respect to my dimension) or do I describe this space efficiently?

Let us be more concrete. In a 2D linear space, I only need two basis vectors to describe all other vectors, by linear combination. Wait, not any tow vectors. Two that are not co-linear. Two (and more vectors) that are colinear, or linearly dependent, have lower rank.
And we do not care about the size of the vectors. So basically, the rank is one mathematical object that both:

check how much a set of vectors is not colinear
does not really bother about their actual size (if $\vec{a}$ and $\vec{b}$ are not colinear, $\vec{a}$ and $10^{-6} \vec{b}$ ae not as well)

So, the issue is not about being big (as scalar factors barely matter), but being "spread" (or dense).

I'll be back latter, building on parts of: Space–Time Codes for High Data Rate Wireless Communication: Performance Criterion and Code Construction,  IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 44, NO. 2, MARCH 1998


Implementing a finite impulse respone (FIR) filter for computing the derivative of a discrete signal

FIR filter with impulse response $h[n] = {-\frac{1}{2}, \frac{1}{2}}$ means that :
$$y[n] = \frac{1}{2}x[n-1] - \frac{1}{2}x[n]$$
This is in some sense a mirror operation of Moving Average of two samples. This is Moving Difference (samples reversed). A High Pass Filter.
So, each output sample is the difference of current input sample with previous input sample. And, this is the result of convolution operation between $x[n]$ and $h[n]$.
But I am not sure that this is will correctly give the Derivative, since for a positively increasing discrete ramp, you would expect the discrete derivative to give positive result. Because for a increasing discrete ramp $x[n]>x[n-1]$. And that FIR filter will give negative result and in addition to negative slope it will also scale the result by $\frac{1}{2}$, which is unwanted.
Correct definition of a Discrete Derivative would be in my understanding, the following:
$$\dot{x} = x[n] - x[n-1]$$
And this can be implemented using FIR filter with impulse response $h[n] = 1,-1$


How to perform 2D image rotation by FFT based image translation using matlab?

The following is the completed code that yields the expected result:

close all; clear all; clc;

I = double(imread('Cameraman.tif'))/255;
figure,imshow(I)

[M, N, s] = size (I);

temp=ones(2*M,2*N);
temp(floor(M/2):floor(M/2)+M-1,floor(N/2):floor(N/2)+N-1)=I;
I=temp; clear temp;

figure,imshow(I)
[M N s]=size(I)
teta = pi/3;
a = tan(teta/2)
b = -sin(teta)

 Nx = ifftshift(-fix(M/2):ceil(M/2)-1);
 Ny = ifftshift(-fix(N/2):ceil(N/2)-1);
 Ix=zeros(length(Nx),length(Ny));
 for k=1:M
   Ix(k,:)=real(ifft(fft(I(k,:)).*exp(-2*i*pi*(k-floor(M/2))*Ny*a/N)));
 end

 Iy = zeros(size(Ix));
 for k=1:M
   Iy(:,k)=real(ifft(fft(Ix(:,k)).*exp(2*i*pi*(k-floor(M/2))*Ny*b/N)'));
 end

 If = zeros(size(Iy));
 for k=1:M
   If(k,:)=real(ifft(fft(Iy(k,:)).*exp(-2*i*pi*(k-floor(M/2))*Ny*a/N)));
 end

 figure,imshow(Ix);
 figure,imshow(Iy);
 figure,imshow(If);

So this works as expected from a 2D rotation about z-axis. The paper argues that this is a fast and accurate method, but I doubt that, as a simple 2D rotation could easily be implemented using raw trigonometric computations and a bilinear interpolation (or bicubic if a sharper result is desired) other than forward and inverse FFTs.

I can say that the method is slow (on matlab, it can be further vectorized however). On the other hand the perceptual quality of the resulting rotated image is fairly good except some little ringing artefacts towards the edges.

I suggest you the ch-5 of Principles of Computer Graphics 2E by FOLEY for further understanding of the mathematics used in these algorithms and applications.


Windowing after filtering a discrete signal to get it back to its original length

You cannot "really" apply a window whose length is different from the size of the signal, since it resorts to a sample-wise product. However, you can for instance apply a composite window, made of a standard window $w$ for the first $N$ samples, and filled with zeros for the remaining samples (to the right). If $w$ is a rectangular window, you get with conv.m what you would have with filter.m, plus the flat tail, as explained by @MattL; the red lines (filter.m) and blue crosses (conv.m) are matching perfectly:

Of course,  you could use a non-rectangular window (here the Hann), and you get a nicer apodization at the edges (the signal is flattened to $0$), but you could as well window the output of filter.m


In the STFT, you apply windowing and Fourier transform on the signal using sliding patches and then combine the resulting transforms, which will help you eventually end up with a uniform time/frequency representation of the signal.

In the wavelet transform case, you apply a filter bank on the overall signal at once. In this way, you obtain a coarse-to fine resolution pattern on the time/frequency representation.

Both methods result in similar time/frequency representations which can be derived from each other.

The major differences:
(1) STFT is uniform yet CWT is not.
(2) You apply STFT on patches, but you apply CWT on the overall signal.
(3) STFT involves Fourier transforms but CWT only requires an orthogonal filter bank.


Calculate symbol rate for QAM modulation

There is no inherent symbol rate for any modulation type.  64 QAM symbols transmit 6 bits per symbol, and 256 QAM symbols transmit 8 bits per symbol.  If the bit rate were to stay constant, then the 256 QAM symbol rate would be 6/8 of the 64 QAM symbol rate, but it appears that that is not the case (assuming your numbers are correct).  That may be because they are using more forward error correction (FEC) to mitigate the higher power requirements of 256 QAM.

Given that the bit rates are not the same, we can't guess what the 128 QAM symbol rate would be.

The general equation for a given modulation type and its symbol rate is:
$symbols/s = \frac{bits/s}{bits/symbol}$

Determining the bits per symbol is usually pretty easy.  It's log2 of the number of possible symbol points.  Thus, the bits per symbol of 256QAM (which has 256 possible points) is 8.

The trick is finding the bits per second.  Wireless transmitters often have the following blocks:

Data source => Compression => Scrambling => FEC Encoding => Pulse Shaping => Modulation

To know the final bit rate you need to know the bit rate of the data source, how much the compression reduces the data rate (if there is compression), and how much the FEC encoding increases the data rate.


Lesser bins histogram looks more equalised than higher bins histogram after histogram equalisation?

I'm not sure if this is the best answer to your question out there, but I would disagree that your 256-bin histogram isn't equalized. Think of it this way.

Each bin represents the count of pixels of a certain intensity. When you're performing histogram equalization, you're moving that bin to a new intensity value, according to the empirical distribution (i.e. histogram of the original image). You're shifting all of those pixels to a new value, which means that if 30% of all your pixels had a certain value, all 30% will again be in the same bin, just with a different value.
Histogram equalization tries to equalize the density of pixel values within any neighborhood. For example, since that huge bin with 30% of all pixels cannot be made smaller (all values of those pixels change to the same new value), equalization will try to move all other bins away from it, so that that neighborhood around the big bin frees up a bit. This is exactly what happens in the picture on the right. Notice that all tall bins have large spacing between them, while all short bins are stuck close together. This ensures that on average you have the same amount of pixels "per unit length of histogram's x-axis" if you will. This brings us to the last point.
When you're using less bins, (10 in your case), you're basically summing or averaging bins of your high-resolution histogram into a low-resolution histogram. Remember that the amount of pixels in a certain intensity neighborhood is what's really equalized, and that really manifests itself when you're summing or averaging components in one of those neighborhoods. Because taller bins have more space between them, and shorter bins have less spacing, you're getting approximately equal height bins when you sum or average. This shows that histogram equalization works as intended.

Hope this helps.


Hint: Expand $e^{j \omega}$ using Euler's formula. This will enable you to express $H(\omega) = R(\omega) + j I(\omega)$ where $R$ and $I$ are the real and imaginary parts and then the magnitude of the frequency response is just $\sqrt{ R(\omega)^2 + I(\omega)^2}$.


Cascaded one-pole sections and relationship to time constant

The difference equation of your single pole filter is

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

where $n$ is the time (sample) index, $x[n]$ is the input, and $y[n]$ is the output. Its transfer function is

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

If you cascade $N$ identical single pole filters, the total transfer function is

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

The inverse $\mathcal{Z}$-transform of $(3)$ is the impulse response of the total system:

$$g[n]=(1-b)^N\frac{(n+1)(n+2)\ldots (n+N-1)}{(N-1)!}b^n,\quad n\ge 0,\quad N>1\tag{4}$$

From $(4)$ it's straightforward to compute the desired value of $b$ for $g[n]$ to decay to $0.01$ of its initial value in a given time $T$. The number of sample values corresponding to time $T$ is $K=\text{round}(Tf_s)$, where $f_s$ is the sampling frequency. We now need to solve

$$g[K]=0.01\cdot g[0]$$

for $b$, which results in

$$b=\left(\frac{0.01\cdot (N-1)!}{(K+1)(K+2)\ldots (K+N-1)}\right)^{1/K}\tag{5}$$

For $K=5$, as in your example, and for $N=2$ cascaded filters you get from $(5)$ a value of $b=0.27821$. For $N=5$ you would need to decrease this value to $b=0.15133$.


The most likely possibility is that you're thinking of partial response signaling, of which GMSK is a special case. This is covered in many textbooks, e.g. in Proakis.


Is it possible to calculate a z-transform for a filter calculated with Parks-McClellan?

If your impulse response is $h(n)$ then the $\mathcal{Z}$-transform of this impulse response (i.e. the transfer function of the filter) is

$$H(z)=\sum_{n=0}^{N-1}h(n)z^{-n}$$

where $N$ is the filter length, i.e. the number of coefficients.

Note that since it is an FIR filter, all its poles are at the origin of the $z$-plane. The zeros are of course the zeros of the polynomial with coefficients $h(n)$.


Audio fingerprinting - 2D (time - frequency) filter / Energy differences

If you split the signal in 33 Frequency bands for every frame you will have a vector band containing all subbands, all right, and now to compute the Energy Difference from each Band you just need to do:

$$E(x)=\sum_{k=1}^{32}band(k) - \sum_{k=1}^{32}band(k+1)$$

For the first frame compute the equation above, and for all others frames block, You can do something like:

$$E2(x)=\sum_{k=1}^{32}band(k) - \sum_{k=1}^{32}band(k+1)$$

$$F(n,m)= 1 if(E2 -E)>0$$
$$F(n,m)= 0 if(E2 -E)<=0$$

$$E(x) = E2(x)$$

This will generate a 32 bits array (0's and 1's) in $F$, Your end result should be a matrix of dimension 32x256, and it is your fingerprint !

EDIT

Answering your other question, if you are following the article exactly as described, the frame rate used in the paper is 5000hz, then each sub-fingerprints 32x1 is equivalent at 0,0116ms of audio its give you an analyze of 58 Frames = 0,0116 * 5000, it does not contain enough information to make a fingerprint, you'll need something comes close to 3 seconds and 256 comes very close look, 0,0116 * 256 = 2,9696, Ok you need do one fingerprint for ever 256 sub-fingerprints and it will create a lot of 32x256 matrix for your audio :-)


Typically, phase ambiguity arises in the demodulation process when the carrier recovery is done via a loop that tracks the fourth harmonic of the received signal.  So why can't the loop track the carrier frequency directly? Well, the received signal doesn't have a carrier signal -- it has a modulated carrier signal whose phase is going to change at every symbol interval which makes the phase untrackable by a typical loop -- and it is necessary to take the fourth power to eliminate the modulation and produce a signal with unchanging phase. As a result, the local oscillator (LO) in the receiver might be in any of four different phases differing by $\pi/2$, and the local reference signals to the mixers in the I and Q branches of the receiver might end up being not the expected
\begin{align}\big(\cos(\omega t), -\sin(\omega t)\big)& & \text{LO and received signal are in phase}\end{align}
as reference signals (note that the Q branch reference leads the I branch reference by $\pi/2$) but instead
\begin{align}\big(\sin(\omega t),\cos(\omega t)\big)& & \text{LO lags received signal by }\pi/2\end{align}
or
\begin{align}\big(-\cos(\omega t), \sin(\omega t)\big)& & \text{LO lags received signal by }\pi\end{align}
or
\begin{align}\big( -\sin(\omega t),-\cos(\omega t)\big)& & \text{LO lags received signal by }3\pi/2.\end{align}
The result is that instead of the I and Q bits appearing at their respective demodulator outputs (the receiver produces $(I,Q)$), the receiver might be producing $(\bar{Q},I)$, or $(\bar{I},Q)$, or $(\bar{Q},I)$ where the overbars denote bit complements.  Differential demodulation is commonly used to eliminate the phase ambiguity.


Fourier Transform of Alternating Periodic Rectangular Pulse

The answer is simple.
I will give 3 points to solve it:

The Fourier transform is linear. Hence $ \mathcal{ F } \left\{ \alpha f \left( x \right) + \beta g \left( x \right) \right\} = \alpha \mathcal{ F } \left\{ f \left( x \right) \right\} + \beta \mathcal{ F } \left\{ g \left( x \right) \right\} $.
Shift in time $ f \left( x - {x}_{0} \right) $ equals multiplication by $ {e}^{-j \omega {x}_{0}} $ in Fourier domain.
Instead of solving for the case above, think of the case you have 2 rectangular signals with twice the period with one multiplied by $ -1 $ and shifted.

In your case, just have a look on the signal of the Positive Pulses. It has a period of $ 4T $.
You have another signal. The signal of the negative pulses.
It is basically the same signal as the positive one (It also has a period of $ 4T $) yet it is shifted by $ 2T $ and as it is multiplied by $ -1 $.
So if you know the transform of the positive one, follow my above points and you have the transform of the negative one and their sum.


What is the advantages of BPSK+OFDM over usual BPSK (with no OFDM Modulation)?

In a pure AWGN channel there is no benefit from OFDM because the noise will affect the received signal in the same way regardless of whether sent as a higher bandwidth serial stream or lower rate parallel streams. Similarly, there isn't much benefit in a flat fading environment either because, by definition, each constituent stream of the OFDM signal will be affected by the same amplitude change and linear phase shift.

The main benefit from OFDM comes with frequency selective fading. That is because the effects of frequency selective fading worsen as signal bandwidth increases. With OFDM, each parallel low rate signal experiences a flat fading channel instead a frequency selective channel at the full bandwidth - i.e. each individual channel can be represented by a single complex value instead a multi-tap filter. That means frequency selective equalization is easily addressed, which allows us to scale to much higher bandwidths - and everybody wants higher bandwidth.

Note that the question refers to BER only, so other OFDM related characteristics such as spectrum mask or resource allocation are left out.


Perfect reconstruction ( Single Channel )

HINT: Note that if

$$Y(z)=X(z^2)$$

then obviously

$$Y(z^{1/2})=X(z)$$

but also

$$Y(-z^{1/2})=X((-z^{1/2})^2)=X(z)$$

from which the desired result follows in a straightforward way.


Accuracy of IFFT(FFT(X)) unit transformation

Your mileage may vary, but in multiple environments on my x86_64 bit linux machine, I get much less error than 1 / 10^9.

julia:

julia> a=randn(32768); f=ifft(fft(a)); print(maximum(abs(f-a)))
1.782802105510316e-15

octave:

octave:1> a=randn(32768,1); f=ifft(fft(a)); disp(max(abs(f-a)))
1.8184e-15

numpy:

In [1]: %pylab
Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

In [2]: a=randn(32768); f=fft.ifft(fft.fft(a)); print(max(abs(f-a)))
7.58800755928e-15

The default for all these environments is float64.  If you want to share how the signal is bandlimited, that might shed more light on the problem.

UPDATE: To satisfy my curiosity, I tried this with a single sine-wav and basically get the same results:

julia> th=collect(0:32767)*(1000*pi/32767); a=sin(th);
       f=ifft(fft(a)); print(maximum(abs(f-a)))
9.393934579042684e-16

Again, without more knowledge about the signal or the computation environment, this is hard to address.


What is 'antenna-chain permutation' and its advantages?

I'll highlight terms that are worth googling in cursive below:

What you're looking at is a diversity receiver.

In this case, the idea is that by combining the signals of multiple antennas, the SNR can be maximized, minimizing the BER.

There's multiple methods of doing so – simply selecting the antenna with the highest RSSI (being probably the one with the best SNR) is called selection combining. When thinking about outage probabilities, this leads to a Bernoulli distribution, ie. the more antennas you have, the less like it is that your communication fails totally (but the gain of having more antennas quickly decreases after the second).

Another method is just adding up the signals – you can do that, and you'd get what's called equal gain combining. You'd get some sort of "oversampling" gain, supressing noise that is independent in each receiver branch, but a strong interferer might distort the result.

Then, you can add the receive signals up, but weigh them by some function that is proportional to their "quality". That way, you get maximum ratio combining, at least if the weights are basically proportional to SNR.

In a real-world receiver with a limited amount of ressources dedicated to assessing the signals, a simple ranking in three tiers of signal quality does sound like a feasible approach. For example, spending much power on synchronizing the strongest signal (i.e. estimating the phase of the channel impulse response) makes sense, since an error here has the gravest effect, whereas simpler approaches for the signals that are used to "push down" uncorrelated noise power (maybe aided by the strongest' chains info) might still make sense. Notice that modern WiFi receivers are pretty mighty beasts – not what you'd need to feel bad about if you don't understand them at once.

Another explanation might be much more in the analog than in the digital domain: For the strongest signal, you might want to use a different receive amplifier than for the weakest one – amplifiers are either a) low-noise, or b) high dynamic range or c) power-efficient (not to mention d) through z) cheap).