If you need to discriminate between 60 and 65 Hz signals, the fact that you're sampling at 25 kHz means that you're basically doing real time signal processing: the Nyquist Frequency is far above the frequencies of interest.

Thus, from classical Fourier theory, you're going to need on the order of 1.0/5 Hz = 0.2 seconds of data, no matter how you process it.  If you do FIR, that means that you'll need at least 5000 coefficients --- and probably more like the 25000 that you mention.

It's possible that you could use much smaller IIR filters.  You'll do less math, but you'll still have to chew through at least 5000 samples for fundamental Fourier reasons.  And you'll need to be careful with that IIR filter, because the pole will be very close to the unit circle: you'll probably want to convert your data to float or double rather than an integer format (which would likely work for FIR).

Now, if there is some prior knowledge you have about your signals --- "high SNR, knowledge of the phase of the 60 and 65 Hz signals..." then you can create discrimination algorithms that are faster.  But all you know is 60/65 Hz, you're sort of out of luck.


When a continuous time signal is sampled in discrete intervals, the effect in the frequency domain is that the CTFT becomes periodic in frequency - that is, you end up with "copies" of the original CTFT, with each copy centered at frequencies equal to k*1/T (k = any integer).

So, when you sample at 12/pi you will get a copy of the original CTFT centered at w=0, w=12/pi, w=-12/pi, w=24/pi, w=-24/pi, and so on.

Since you have frequency content between 2<|w|<4 in the w=0 case, the w=k*12/pi cases will have frequency content between 2<|w-k*12/pi|<4.

Be careful - if frequency content from multiple copies have frequency content in the same places on the frequency axis, the effect is that they add together.  This is called "aliasing" and is usually an unwanted effect of sampling.

Here is a nice graph that shows the effect of sampling a continuous time signal on the FT.


Hidden Markov Models - Distinct Observation Symbols

Firstly, this is not because you have two words to identify that you need $N = 2$ states. Your goal is not to train a model with two states - one for each word to recognize - but to train 2 models, one for each word to recognize, and each of these models will have as many states as necessary. In fact, each state in your HMM should correspond to a distinct "stage" in the pronunciation of a word - and will very likely correspond to a phoneme. Your vocabulary size (here, two: "stop" and "go") is external to this. For "stop", there are 4 phonemes. For "go", there are 2 phonemes. So you train a 4-state left-to-right model on the "stop" data; and independently of this, a 2-state left-to-right model on the "go" data. To recognize a word given its MFCCs, you evaluate which of these two models has the highest likelihood given the data. If you had to recognize words within a lexicon of 10 words, you would similarly train 10 HMMs, one for each word, each of these models having a number of states suitable to the length/complexity of the word to recognize.

You need to step back and ask yourself "why HMMs in the first place?". We need HMM for speech recognition because words are made of a sequence of distinct elements in sequence (phonemes). If we want to describe/recognize the word "stop", we need to learn a description which is expressive enough to capture that "first it sounds like ssss, for a short while, then it is tttt for a short while, then it is oooo for a longer amount of time, then it is pppp for a short moment". HMMs are a good match for expressing that - states are phonemes, the transition matrix (which will be here diagonal + upper diagonal) indicates that we move through the word from first phoneme to last phoneme, staying a variable amount of time in each phoneme, and the distribution associated with each state indicates how each phoneme translates into your acoustic features.

It seems also that you are mixing up discrete HMM (in which the observations are drawn from a discrete distribution associated with each state) with continuous HMMs (in which the observations are scalars or vectors, characterized by a continuous distribution such as a gaussian). So the parameter $M$, number of distinct observation symbols, is irrelevant in your case, since your observations are 13-dimensional vectors, an uncountable set! ($M$ would be... the cardinality of the continuum).

I am afraid the introduction material you have picked is not directly relevant to speech recognition - though it is useful for applications in which HMMs are used to recover hidden structure from a discrete observations (and there are many of them, for example parsing/tagging in NLP). Try to master this material without thinking much about your speech recognition problem, and then move on to material about continuous HMMs with multivariate normal distributions - and finally to continuous HMM with mixtures of multivariate normal distributions (since this is what is likely to work best for speech).


Determining pulse shape for Viterbi Equalization

I'll answer my own question with what I've learned in the past 28 days. Hopefully I've got it right!

Viterbi equalization is just one way to equalize the channel, as an alternative to some other method such as convolving the signal with an equalizing filter. Estimating the channel is an entirely separate problem which must be solved independently of the equalization method.

The trouble with an equalizing filter is it compensates for dips in the channel frequency response by amplifying those frequencies. Unavoidably, the noise is also amplified.

There's a tradeoff to be made: the flatter the frequency response is made, the lower the error due to ISI, but the higher the error due to noise. There's an optimal amount of equalization which minimizes the error due to noise and ISI. An equalizer that tries to hit this sweet spot is a minimum mean square error (MMSE) equalizer.

Viterbi equalization offers an alternative solution. The channel estimate will produce an estimated channel impulse response. The coefficients in this impulse response can be used to build the trellis, and the Viterbi algorithm can then find the most likely sequence of transmitted bits, taking into account the estimated distortion from the channel. The receive filter remains a simple matched filter, with no equalization. The Viterbi equalizer is not linear (unlike the equalizing filter alternative), so it's able to compensate for deep spectral nulls in the channel response in a way the equalizing filter can't. The cost is increased complexity.


The purists will tell you it should be called a Hann window, or VonHann, but not Hanning.

From a math perspective, the Hann window is also equivalent to a sine squared window at half the frequency, so by taking the square root you are getting a sine window.

Why use a window at all?  The seminal paper on the topic is:

harris, fred. (Jan 1978). "On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform"

A link to the paper can be found in reference [10] of the Wikipedia article on windows.

Hope this helps.

Ced


How to reduce latency under mean filter for high noise?

Although this is an old entry and probably the author of the question has already solved the issue, I leave here my solution, just in case, it could help someone else.
As explained in the question, you need a kind of real-time "smoother" for the trajectory. You are going to have trade-offs: speed, latency, energy efficiency, etc.
I propose to use Alpha-Beta filters, Moving Average Filters or some kind of Adaptative Average filter. They should be faster, easy to implement and with less latency than Kalman filter or similar. I have made a couple of tests, trying to reproduce the type of data used in the question, and testing those filters with quick(and nasty) own implementations. The results are shown below in the next figures:

Cheers.
P.D. Probably, the vibration noise commented in the original question is not absolutely uncorrelated. The data I have generated is using random noise.


Is there a physical reason why phase should be uniformly distributed?

A zero mean complex white noise process has a uniform phase distribution.  The uniform phase distribution would occur for any complex noise process whose I and Q (real and imaginary) components are zero-mean, independent and identically distributed.
The thermal noise floor common to radio receivers for receiving signals without other specific interference is typically well approximated by a zero-mean complex white noise Gaussian distribution; meaning a white noise Gaussian distribution on I (the real component) and Q (the imaginary component) and these two are independent from each other. This makes sense if you consider that the real and imaginary values for any given sample have no relation to each other and are also therefore independent. Within this restriction of having an independent and identically distributed random process for I and Q, the actual phase of any sample is equally likely as any other phase and is therefore uniformly distributed (regardless if the system is white noise or band-limited).
Complex quantization noise (using two real A/D converters in practice; one for the I channel and the other for the Q channel) will also be well approximated by a white noise process where in this case the phase is uniformly distributed as above but unlike the thermal noise example given above where the signed amplitude was approximated as Gaussian distributed, the amplitude for quantization noise is also uniformly distributed over a quantization level.


A general conceptual understanding of moving average

That system you're referring to is a first-order recursive filter

$$y[n]=\alpha x[n]+(1-\alpha)y[n-1],\qquad 0<\alpha<1\tag{1}$$

where $n$ is the (integer) time index, $y[n]$ is the output, and $x[n]$ is the input. You chose $\alpha=\frac12$.

That filter is also called exponentially weighted moving average. It has theoretically infinite memory and weighs the past input values with exponentially decreasing weights.


You might have to make your own from a 3d face database, which may not be easy to do. This website has a good listing of public databases. http://www.face-rec.org/databases/


Why does a simple delay system with $h[n]=\delta[n-n_0]$ have a linear phase?

In response to "1) What is the filter order of the simple delay system?"

The order of a filter is the power of the highest nonzero coefficient of the Z-transform. For the simple delay system given, the Z-transform is given as $H(z) = z^{-n_0}$ so it has order $n_0$.

In response to "2) Does it have symmetric filter coefficients?"

The coefficients are not symmetric if you stop at $n_0$, but adding additional zero coefficients to the signal will not alter it and by doing so, one may describe an equivalent filter that has coefficient symmetry. For example, if the delay is 2 and our original filter is [0, 0, 1], then we can add two more zeros to obtain [0, 0, 1, 0, 0]. These two filters have identical behavior, but the second one can be said to have symmetric coefficients. I would argue that a simple delay does have symmetric coefficients since it can always be put into an equivalent form that has symmetric coefficients.

In response to "3) Why does it have a linear phase?"

The simplest way to understand why a delay results in a linear phase is to examine the output for a sinusoidal input that is delayed by a constant amount. If we have some sinusoidal signal such as $$x[n] = \cos\left(2 \pi f_0 n\right)$$ and we consider what happens after it is passed through our system $$y[n] = x[n] * h[n]$$ which in our case is simply
$$y[n] = x[n-n_0] $$
or if we substitute in $x[n]$ we obtain
$$y[n] = \cos\left(2 \pi f_0 n - 2 \pi f_0 n_0\right).$$
Notice that the resulting output of the system is a sinusoidal signal with the same frequency as the input, but with an additional phase shift. In this case, the phase shift is $2 \pi f_0 n_0$. However, if we are to consider the phase as a function of the input frequency of the sine wave (that is instead of fixing the frequency as $f_0$, we allow it to vary and rename it $f$), then we obtain a phase function
$$ \phi(f) = 2 \pi n_0 f $$
and from this, you can see that the phase response varies linearly in $f$ with a slope of $2 \pi n_0$ radians. We can use this phase function to predict the resulting phase offset at various frequencies. And this prediction is related to the frequency through a simple linear relationship. That is why we say that it has linear phase.


Detecting multi-path propagation when one acoustic signal is dominant

It seems to me that, in the conditions you describe, the most you will see is flat fading. One indication that flat fading is occuring is that a change in the positions of transmitter and receiver, and/or the reflectors causing the multipath, produce a significant change in the received power. Usually the change has to be of the order of $\lambda/2$ meters, so this may be unfeasible in your case.


Find Amplitude of each period of time series

I assume from your solution that your time series closely fits a sine wave.  The method you are using is quite good if you restrict yourself to the set of data which is near the peaks.

The more standard way is to find a best fit sinusoid using linear algebra techniques.  You have stated the period, so you know the frequency.  Let's assume you don't know where the zero crossings are.  The problem then becomes to find the best values of $(a,b)$ so that

$$ Y = a C + b S $$

Where $Y$ is your signal values as a vector, $C$ is a cosine curve, and $S$ is a sine curve over one period.  Now dot this equation with $C$ and then $S$ to get:

$$ Y \cdot C = a C \cdot C + b S \cdot C $$
$$ Y \cdot S = a C \cdot S + b S \cdot S $$

The dot products are scalars.  Since $S$ and $C$ are orthogonal, their dot product is zero.  The dot product of $S$ or $C$ with itself is $N/2$ where $N$ is the sample count.  The equations then become

$$ Y \cdot C = a \cdot N/2 $$
$$ Y \cdot S = b \cdot N/2 $$

Solve for $a$ and $b$.

$$ a = \frac{ Y \cdot C }{ N/2 } $$
$$ b = \frac{ Y \cdot S }{ N/2 } $$

Your amplitude can now be found from $a$ and $b$.

$$ A = \sqrt{a^2+b^2} $$

This method is equivalent to a single bin of a Discrete Fourier Transform (DFT), so if you follow this then you understand how a DFT works.

To calculate the dot product, you take a summation:

$$ Y \cdot C = \sum_{n=0}^{N-1} { Y[n] C[n] } $$
$$ Y \cdot S = \sum_{n=0}^{N-1} { Y[n] S[n] } $$

That's all there is to it.


Why are numerator coefficients named $b$ and denominator coefficients named $a$?

My suspicion is that this ordering comes from the difference equation, which in most texts precedes the $\mathcal Z$-transform:

$$a_0 y[n] + a_1 y[n-1] + \cdots = b_0 x[n] + b_1 x[n-1] + \cdots$$


What do you call it when you add zeros between the elements of a vector?

In the context of filter banks, simply upsampling (without filtering, or with the trivial all-pass filter, as Wikipedia relates it to interpolation), denoted by an up-arrow, as illustrated below fro a 2-fold or order-2 upsampling:

This is described in Filter banks: Decimation, Interpolation, and Modulation, and the  interpolation filter associated with upsampling is called synthesis filter.

As noted by @Jason R, one also finds "expansion". Expansion is also found as a streching for continuous signals, while compression is the effect in the dual Fourier domain. Sometimes, one can find "zero-insertion" or "zero-stuffing".

Actually in image processing some use "upsampling" as a proxy for pixel increase with filtering.


What does it mean by DN in salt-pepper noise?

From the same course, you can find the following indication in an earlier set of slides:

  Input quantized image pixel values (integers): Digital Number (DN)

From this you can infer that

  Salt noise: DN = maximum possible

indicates that salt noise shows up as pixels with the largest possible pixel values (e.g. 255 if the picture is a grey scale picture represented by pixel in the [0-255] range), and that

  Pepper noise: DN = minimum possible

indicates that pepper noise shows up as pixels with the smallest possible pixel values (e.g. 0 if the picture is a grey scale picture represented by pixel in the [0-255] range).


How do I efficiently evaluate a convolution integral between different-sized domains?

Essentially, to my understanding this problems boils down to a fast-convolution problem:

You have a function $g(x)$, which has a small domain and you have a function $K(\tau)$, which has a big domain and you want to calculate

$$f(x) = \int_{-\infty}^{\infty}g(x')K(x-x')dx'$$

in a discretized fashion, where you want to calculate as efficient as possible. Especially, you are interested in $f(-L/2+nL/N)$ for $n=0,\dots,N-1$.

When considering the discretization you need to take care of the bandwidth of $g$ and $K$ and you need to sample with the according Nyquist frequency at least. Let's assume that $g$ can be sufficiently sampled with $g(-\frac{a}{2}+n\frac{a}{N})$ (i.e. the sampling is fast enough regarding Nyquist frequency). Let's also assume, that $K(\tau)$ has the same or smaller bandwith, i.e. the sampling rate for $g$ is also sufficient for $K$.

Then, let's further assume $L=qa, q\in\mathbb{N}$, i.e. L is a multiple of $a$. Define the following functions:

$$g[n]=g(-\frac{L}{2}+n\frac{L}{qN}), \quad n=0,\dots,qN-1$$
and
$$K[n]=K(-\frac{L}{2}+n\frac{L}{qN}), \quad n=0,\dots,qN-1$$

Note, that $g[n]$ is zero for most $n$, but that is fine here.

Now, you have

$$f[n]=\sum_n'g[n']*K[n-n'], \quad n=0,\dots,qN-1$$

Which is a higher-resolution version of your original $f_n$ (It has qN instead of only N samples within the interval $\pm\frac{L}{2}$). So, what you actually want to have is a downsampled version of this $f[n]$ given by

$f_n=f[qn]$ (where $f_n$ is the $f_n$ from your post, i.e. with $n=0,\dots,N-1$.

So, this is also a discrete convolution problem, and it can be solved via FFT (note that the FFT actually considers circular convolution, if you don't zero-pad). However, since your $g[n]$ has a much smaller domain that $K[n]$, you can resort to fast-convolution algorithms (Overlap-save or Overlap-add).

To summarize: The trick is actually just to sample both signals with the same sampling frequency (accepting that $g[n]$ contains a lot of zeros) and then understanding this problem as a simple discrete convolution.


You will get

$Y(z) = \frac1M\sum_{m=0}^{M-1} X(e^\frac{-j2\pi{m}}M*z)$

If you were to upsample first and then downsample you would get

$Y(z) = \frac1M\sum_{m=0}^{M-1} X(e^{-j2\pi{m}}*z)$

which just simplifies to (since $e^{-j2\pi{m}}=1$ for $m\in\mathbb{Z}$)

$\frac1M\sum_{m=0}^{M-1} X(e^{-j2\pi{m}}*z) = \frac1M\sum_{m=0}^{M-1} X(z) = X(z)$

Basically, if you upsample and downsample you can perfectly reconstruct the signal but if you downsample and then upsample, you are aliasing in the summation over rotations in the z-domain of the input.


CDMA : Calculating M from m1 * w1, m2*w2, etc M as Walsh Code, M as Orthogonal COde

You could have given a bit more background and explanation in your question, but as far as I understand your question, $m2$ is your data and $w2$ is the spreading code. The spreading ratio is 4 because the length of the spreading sequence is four times the data length. So you have to stretch out your data by a factor of 4 to line it up with the spreading sequence. Then you can do the multiplication as follows:

$$\begin{array}{rrrrrrrrrrrrrrrrr}m2: & 1 & 1 & 1 & 1
& -1 & -1 & -1 & -1 & -1 & -1 & -1 & -1 & 1 & 1 & 1 & 1\\
w2: & -1 & -1 & 1 & 1& -1 & -1 & 1 & 1& -1 & -1 & 1 & 1& -1 & -1 & 1 & 1\\
\text{result}:
& -1 & -1 &1&1&1&1&-1&-1&1&1&-1&-1&-1&-1&1&1\end{array}$$


Generating good synchronization sequences

A synchronization sequence generally needs the property that its autocorrelation function resembles an impulse.  There are two possible autocorrelation functions that can be considered. For a (real-valued) sequence $x$ of length $N$, the periodic autocorrelation function is
$$R_x[n] = \sum_{k=0}^{N-1}x[k]x[k+n]$$
where the sequence is assumed to extend periodically so that
$x[k+n] = x[k+n\bmod N]$.  In this case, one commonly used solution to
choose $x$ to be a maximal-length linear binary feedback shift-register
(LFSR) sequence of
period $N = 2^m-1$ whose periodic autocorrelation function is
$$R_x[n] = \begin{cases} N, & n \equiv 0 \bmod N,\\ -1, & n \not\equiv 0 \bmod N,
\end{cases}$$
which resembles an impulse train of period $N$ except that the "out-of-phase"
value of the periodic autocorrelation function is not identically
zero, but is small (and constant)
compared to the "in-phase" value of $N$.

The other form of autocorrelation function is the aperiodic autocorrelation
function which is given by
$$C_x[n] = \begin{cases}\sum_{k=0}^{N-1-n}x[k]x[k+n], & 0 \leq n < N,\\
\sum_{k=0}^{N-1+n}x[k-n]x[k], & -N < n < 0,\\
0, & |n| \geq N.
\end{cases}$$
Now, for binary sequences where $x$ takes on values in $\{+1,-1\}$,
the value of $C_x[n]$ alternates between even and odd integers as
$n$ varies from $-(N-1)$ to $(N-1)$, and so a "flat" out-of-phase
aperiodic autocorrelation function is impossible to achieve.
The closest we can hope to come is to find sequences for which
$|C_x[n]| \leq 1$ for $1\leq n < N$.  Sequences with this property
are called Barker sequences, and the longest known Barker
sequence is of length $13$.  it is widely believed, though not
proven, that Barker sequences of longer lengths do not exist,
and numerical studies show that if such a sequence exists, it
it must be so long that it will not be be usable in any practical
scheme as a synchronization preamble.


As pointed out in other answers, the Central Limit Theorem is one reason why Gaussian noise is so important as a model. One other reason that hasn't been explicitly mentioned and that I would like to point out is the fact that Gaussian noise is completely described by second order statistics, which are relatively easy to measure. E.g., if you have a univariate distribution and if it is Gaussian, then you know everything you can ever know about it by knowing (measuring) its mean (first order) and its variance (second order). There are no higher order statistics, which is good because they are much harder to measure reliably. This fact is of course no justification for using a Gaussian model, but it is a very handy property if we think that a Gaussian model is justified.


Hilbert transform in the frequency domain

The discrete-time Fourier transform (DTFT) is always periodic. This is also the case for the frequency response of the discrete-time Hilbert transformer. For this reason, the ideal frequency response is not only zero at DC but also at Nyquist, which corresponds to index $2$ for a signal of length $4$. Consequently, the correct way to do what you're trying to do is:

x = [1,2,7,3];
X = fft(x);
hil = [0,-1i,0,1i];
Y = hil .* X;
y = ifft(Y);
y =

  0.5  -3  -0.5  3


Spectral baseline distortion simulation - Bruker smile

My guess is that somewhere in the processing you are using an AR (or ARMA) model that is being over / under fitted and that is what is causing the issue.

If I make a fake signal like this in R:

group_delay <- 84
T <- 350

omega <- 2*pi*0.252352890
tau <- 0.005
phi <- 2*pi*0.0989038
t <- 1:T

bruker_signal <- rep(0,T)
t_index <- seq(group_delay, T)
bruker_signal[t_index] <- exp(-(t_index-group_delay)*tau)*(
    sin(omega*t_index) + sin(1.5*omega*t_index + phi)) +
     rnorm(length(t_index),0, 0.1)

i.e. it's a noisy, delayed, damped sum of two sinusoids and then do AR modelling of the spectrum with different orders, then I get something like the smiles / frowns you are seeing.

An earlier attempt I made used a low-pass filtered version of the signal added to the same thing, modulated to $f_s/2$.

smile_filter <- butter(3,0.1)
smile_signal <- filter(smile_filter, bruker_signal)
plot(abs(fft(bruker_signal+2*smile_signal + 2*smile_signal*rep(c(1,-1),     lenght.out=T)))[1:175], type='l')
title('Attempt at Smile Spectrum ')

but that seemed way noisier than what you're after.


Does that mean awgn function behavior depends on the zero padding (length of the signal)?

Yes, it does depend on the length of the signal.

How does the 'awgn(input_signal,snr_db,'measured')' function calculates the input signal power?

For an input signal of length $N$, $x[n]$, the awgn function calculates the power as $\frac{1}{N} \sum_{n=1}^N |x[n]|^2$.

How to make the power calculation independent of the length of the input signal?

The above definition defines $N$ as the length of the signal, including the zero valued samples. If you don't want to count those towards your signal length, then you needs to have an effective signal length $N_0$ that only counts non-zero samples.
sigPower = mean(abs(qam_mod_u(qam_mod_u ~= 0)).^2);
noisePower = sigPower/10^(SNRdB/10);
noise = sqrt(noisePower/2)*(randn(size(qam_mod_u)) + 1j*randn(size(qam_mod_u)));
ch_out_interplt = qam_mod_u + noise;

Is there any better way to simulate the awgn noise in MATLAB?

The awgn function is fine as long as you understand what it is doing. For situations when the signal is turning off and on, it might be easier to do the noise on your own for exactly the reason of your question. You can also make it easier by putting this all in your own function too.


Inverse discrete time Fourier transform with differentiation

Hint:
differentiate from Fourier Transform equation:
$$X(e^{j\omega})= \sum\limits_{n=- \infty }^{\infty}x[n]e^{-jn\omega} \Rightarrow \frac{d}{d\omega}X(e^{j\omega})= \frac{d}{d\omega}\sum\limits_{n=- \infty }^{\infty}x[n]e^{-jn\omega}$$
and then find what you want.


You did not implement a recursive filter, but a non-recursive filter. Note that in your implementation each output only depends on possibly delayed input values, whereas in a recursive filter, the output depends on delayed input and output values.

If the coefficients of one second-order section are given by
$[b_0, b_1, b_2, 1, a_1, a_2]$, then if $x[n]$ is the input sequence, the output sequence $y[n]$ is computed as

$$y[n]=b_0x[n]+b_1x[n-1]+b_2x[n-2]-a_1y[n-1]-a_2y[n-2]$$


Sphere decoding for non-binary modulation

This channel model is known as a multiplicative channel, and it differs from the AWGN channel in that the signal amplitude is scaled by a gain factor.

If the gain factor is real (that is, if the channel simply attenuates the signal), then we can write $\Re (y) = G\Re(a) + \Re(w)$ and $\Im(y) = G\Im(a)+\Im(w)$. In other words, the real and imaginary parts of the received signal are orthogonal and they can be detected separately.

In a case like this, using SD is overkill. For instance, if the transmitted constellation is QPSK, the optimum detection problem involves calculation of four distances in $\mathbb{R}^2$, which is not a hard problem. Furthermore, the symmetry of QPSK actually reduces the problem to that of finding in which quadrant $y$ lives in, which involves two sign checks.

To make things more interesting, consider a flat-fading wireless channel, where the gain $G$ is complex. In this channel, $\Re(y)$ depends on both the real and the imaginary parts of $a$, so the problem cannot be neatly separated in two anymore. If we write $G=|G|e^{j2\pi\angle{G}}$, then we can see that the wireless channel rotates and scales the transmitted constellation.

Even more interesting is the case where we have receive diversity, or mulitple receiver antennas operating concurrently. In this case $G$ becomes a vector, each element of which is the complex channel gain between the transmitter and one receiver. Now, $y$ is also a vector. For a large number of antennas and/or a large constellation, the problem of estimating $y$ becomes more and more complex.

Consider the case where we use 32-QAM and 4 receiver antennas. The ML calculation now involves calculating $32^4\approx10^6$ distances! And that is to estimate just a single symbol. It is in this situation that the SD algorithm starts to shine, since it can dramatically reduce the number of combinations to check.

Unfortunately, in a case like this the SD problem cannot be neatly separated into two. However, there are (more complex) ways to break the calculations into independent blocks that can be computed in parallel. Search for "parallel sphere decoder" in Google scholar or similar engines: there are many results that can help you speed up the SD algorithm, and still find an optimum solution (or close to).


You need to use a bank of bandpass filters before applying the Hilbert transform to get per frequency band amplitude and phase change over time.


Gaussian Equivalent of Convolving an Image 50 Times with a Box Filter

What you should do is work with the Variance.
The Variance of Discrete Random Variables with support of $ \left[ -l, l \right] $ is given by (Notice that since the Expected Value is $ 0 $, the Variance is given by the 2nd moment)):
$$ \sum_{k = -l}^{l} \mathbb{P} \left( k \right) {k}^{2} = \sum_{k = -l}^{l} \frac{1}{2l + 1} {k}^{2} = \frac{1}{2l + 1} \sum_{k = -l}^{l} {k}^{2} $$
Convolution of a function which is the PDF of a random variable is equivalent of addition of the Random Variables (Assuming they are independent).
The Variance is linear with addition of Independent Variables hence the Variance of $ 50 $ Convolutions (Which is addition of 51 variables) is given by: $$ \frac{50}{2l + 1} \sum_{k = -l}^{l} {k}^{2} $$
Now generate a Gaussian Blur with $$ \sigma = \sqrt{ \frac{50}{2l + 1} \sum_{k = -l}^{l} {k}^{2} } $$ apply it on the image and it will be the best approximation for applying 50 times the Box Filter with Radius of $ l $.
By the Central Limit Theorem we can get better and better approximation as the number of Kernels gets higher.
Remarks

Usually Gaussian Blur with large STD is approximated by few (4-5) iterations of Box Blur. Yet in this case since there is a large number of convolutions the proper way is to calculate the STD of equivalent Gaussian Blur and then apply it either directly or by accelerated method (For example, approximation by few box blur iterations).
Box Blur can be applied in O(n) manner. Namely complexity is independent of the box kernel radius. It can be done by Integral Image or Running Sum.
Pay attention that if we do 50 Convolution of Box Blur with itself it is equivalent of 51 additions of Random Variables. Yet above the first convolution is done between the Box Blur Kernel and the images hence 50 Convolution means addition of 50 Random Variables.

MATALB Code
boxRadius   = 4;
boxGrid     = [-boxRadius:boxRadius].';
boxSupport  = (2 * boxRadius) + 1;
vB          = ones([boxSupport, 1]) / boxSupport;

boxStd = sqrt(sum((boxGrid .^ 2) .* vB));

vT = vB;

numConv = 9;

for ii = 1:numConv
    vT = conv2(vB, vT, 'full');
end

kernelRadius    = floor(size(vT, 1) / 2);
vKernelGrid     = [-kernelRadius:kernelRadius].';
kernelSupport   = (2 * kernelRadius) + 1;

kernelStd = sqrt(sum((vKernelGrid .^ 2) .* vT));

kernelStd / boxStd %<! Should be 'numConv + 1

vGaussianKernel = (1 / sqrt(2 * pi * kernelStd * kernelStd)) * exp(-(vKernelGrid .^ 2) / (2 * kernelStd * kernelStd));

figure();
plot(vKernelGrid, [vT, vGaussianKernel]);
xlabel('Kernel Support');
ylabel('Kernel Value');
legend({['Convolved Box Blur'], ['Gaussian Blur']});

Closed Form Formula for the Variance of Centered Uniform Discrete Random Variable
As seen above if $ x $ is a centered uniform Random Variable on the support $ \left\{ -l, -l + 1, \ldots, 0, 1, \ldots, l - 1, l \right\} $ then its variance is given by the second moment (As it is centered which means the mean value is $ 0 $):
$$ \operatorname{Var} \left( x \right) = \sum_{k = -l}^{l} \mathbb{P} \left( x = k \right) {k}^{2} = \frac{2}{2 l + 1} \sum_{k = 1}^{l} {k}^{2} $$
The above sum of squares can be solved either by integral or nice trick form the Binomial Identity:
$$ {\left( k - 1 \right)}^{3} = {k}^{3} - 3 {k}^{2} + 3 k - 1 \Rightarrow {k}^{3} - {\left( k - 1 \right)}^{3} = 3 {k}^{2} - 3 k + 1 $$
Summing each side from $ k = 1 $ until $ n $:
$$ \sum_{k = 1}^{n} \left( {k}^{3} - {\left( k - 1 \right)}^{3} \right) = {n}^{3} = 3 \sum_{k = 1}^{n} {k}^{2} - 3 \sum_{k = 1}^{n} k + \sum_{k = 1}^{n} 1 = \sum_{k = 1}^{n} {k}^{2} - 3 \frac{n \left( n + 1 \right)}{2} - n $$
Rearranging yields:
$$ \sum_{k = 1}^{n} = \frac{n \left( n + 1 \right) \left( 2 n + 1 \right)}{6} $$
Hence, for $ x $ as above:
$$ \operatorname{Var} \left( x \right) = \sum_{k = -l}^{l} \mathbb{P} \left( x = k \right) {k}^{2} = \frac{2}{2 l + 1} \sum_{k = 1}^{l} {k}^{2} = \frac{2 l \left( l + 1 \right) \left( 2 l + 1 \right)}{6 \left( 2 l + 1 \right)} = \frac{l \left( l + 1 \right)}{3} $$


Why is peak detection in chromatography not completely automatic?

This is a problem we are also currently working on (added two recently published papers below), because we found it unresolved yet, and high-throughput experiments nowadays require more automated high-throughput analyses.
A first set of causes are DSP-like disturbances on the "interesting peak data": non-linearity (like saturation), added noise (with unknown distributions) and baselines (with unknown regularity). They can related to solute or matrix effects, type of chromatography (liquid, gas), columns, detectors, the complexity of the mixture, etc. We can include more or less visible processing performed inside the measurement tools: sampling (how any points per peaks?), smoothing, normalizations.
In gas chromatography there are models like Gaussian shapes, that can be explained by physical models (plates). However, separated peaks can exhibit distorsions such as fronting and tailing, or other asymmetries, that can vary throughout the chromatogram in uneven manners.
Even with "nicely"-sampled discrete Gaussians, very close peaks may be difficult to separate or distinguish (shouldering effects). More globally, deconvolution of sums of Gaussians remains difficult, because the direct problem is ill-posed: there are many ways to fit a single Gaussian with one of several Gaussians, and algorithms require additional optimization terms (constraints, regularization), that come with tuning parameters when combined with the data fidelity term (often a least-squares, but other choices are interesting).
As far as I know, such issues exists in other analytical chemistry modalities: NMR, infrared, mass spectrometry, Raman etc., each with their own issues. Variable peak shifting related to molecule interactions in NMR is a challenge for subsequent (blind) source separation algorithms.
There are works on using more machine learning techniques, yet it is likely that using physical or DSP models would help them.
A possible starting point on the difficulties of fitting: Dependence of Chromatogram Peak Areas Obtained by Curve-Fitting on the Choice of Peak Shape Function, 1997
We recently considered a couple of methods towards reconstruction or recovery of physico-chemical signals. They definitely are not the end of the game:
SPOQ ℓp-Over-ℓq Regularization for Sparse Signal Recovery applied to Mass Spectrometry, IEEE Transactions on Signal Processing, 2020 (Matlab code SPOQ toolbox)

Underdetermined or ill-posed inverse problems require additional
information for \ldd{d} sound solutions with tractable optimization
algorithms. Sparsity yields consequent heuristics to that matter, with
numerous applications in signal restoration, image recovery, or
machine learning. Since the ℓ0 count measure is barely tractable, many
statistical or learning approaches have invested in computable
proxies, such as the ℓ1 norm. However, the latter does not exhibit the
desirable property of scale invariance for sparse data. Extending the
SOOT Euclidean/Taxicab ℓ1-over-ℓ2 norm-ratio initially introduced for
blind deconvolution, we propose SPOQ, a family of smoothed
(approximately) scale-invariant penalty functions. It consists of a
Lipschitz-differentiable surrogate for ℓp-over-ℓq quasi-norm/norm
ratios with p∈]0,2[ and q≥2. This surrogate is embedded into a novel
majorize-minimize trust-region approach, generalizing the variable
metric forward-backward algorithm. For naturally sparse
mass-spectrometry signals, we show that SPOQ significantly outperforms
ℓ0, ℓ1, Cauchy, Welsch, SCAD and Celo penalties on several performance
measures. Guidelines on SPOQ hyperparameters tuning are also provided,
suggesting simple data-driven choices.

Sparse Signal Reconstruction for Nonlinear Models via Piecewise Rational Optimization, Signal Processing, 2021

We propose a method to reconstruct sparse signals degraded by a
nonlinear distortion and acquired at a limited sampling rate. Our
method formulates the reconstruction problem as a nonconvex
minimization of the sum of a data fitting term and a penalization
term. In contrast with most previous works which settle for
approximated local solutions, we seek for a global solution to the
obtained challenging nonconvex problem. Our global approach relies on
the so-called Lasserre relaxation of polynomial optimization. We here
specifically include in our approach the case of piecewise rational
functions, which makes it possible to address a wide class of
nonconvex exact and continuous relaxations of the ℓ0 penalization
function. Additionally, we study the complexity of the optimization
problem. It is shown how to use the structure of the problem to
lighten the computational burden efficiently. Finally, numerical
simulations illustrate the benefits of our method in terms of both
global optimality and signal reconstruction.


Is the Laplace transform a special case of Fourier transform? (Not the other way around)

The Fourier Transform is the Laplace Transform with the complex variable s restricted to be the imaginary axis on the s plane. For this reason the Fourier Transform only exists when the imaginary axis is within the region of convergence. The variable s is called a "complex frequency" as it is the frequency variable that can take on real ($\sigma$) and imaginary ($\omega$) components. That said, I would view the Fourier Transform as a subset of the Laplace Transform, or the Laplace Transform as an expansion on the Fourier Transform that provides a lot more functionality and can exist when the Fourier Transform can't.

This is also the reason that the frequency response for a system with a general transfer function $H(s)$ is given as $H(j\omega)$.

When a system is restricted to $s= j\omega$ as the input, then the input is restricted to be only sinusoids or signals given by $e^{st}$ with $s = j\omega$ which maintain a constant magnitude with time. By allowing s to have real and imaginary components as in $s = \sigma + j\omega$ then we also allow the input to grow or shrink with time, depending on which point in the s plane is used as the input to the system.


As mentioned the Haar filters use a lookup table with the say $A_{i,j} = \sum_{k=0}^{i} \sum_{l=1}^{j} I_{k,l}$
Now consider a the convolution with a filter $H \in \mathbb{R}^{m \times n}$, such that $H_{i,j}=1$, the convolution will can be calculated as
$$\begin{eqnarray}Y_{i,j} &=& \sum_{k=1}^{m} \sum_{l=1}^{n} H_{k,l} I_{i-k+1, j-k+1} \\
&=&  \sum_{k=1}^{m} \sum_{l=1}^{n} I_{i-k+1, j-k+1} \\
&=&  \sum_{k=i-m+1}^{i} \sum_{l=j-n+1}^{j} I_{k, l} \\
&=&  \sum_{k=i-m+1}^{i} \left( \sum_{l=0}^{j} I_{k, l} -  \sum_{l=0}^{j-n} I_{k, l}\right) \\
&=&  \sum_{k=0}^{i} \left( \sum_{l=0}^{j} I_{k, l} -  \sum_{l=0}^{j-n} I_{k, l}\right) - \sum_{k=0}^{i-m+1} \left( \sum_{l=0}^{j} I_{k, l} -  \sum_{l=0}^{j-n} I_{k, l}\right)\\
&=&  \left(A_{i,j} -  A_{i,j-n}\right) - \left( A_{i-m,j} -  A_{i-m, j-n}\right)\\
\end{eqnarray}
$$
Of course you can generalize this for $H_{i,j} = c$ for constant $c$ (by multiplying the resulting image), and you can use the shifting properties of the convolution if the filter is padded.
Now looking to your $H$, it can be decomposed in two rectangular regions constant regions. It can be decomposed as
$$H_1 + H_2 = \begin{bmatrix} 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 & 0\end{bmatrix} + (-1) \begin{bmatrix} 0 & 1 & 1 & 0 & 0 \\  0 & 1 & 1 & 0 & 0  \\  0 & 1 & 1 & 0 & 0 \end{bmatrix}$$
To convolve with a the non-zero submatrix of $H_1$, a matrix in $\mathbb{R}^{3 \times 1}$ you will use $A_{i,j} + A_{i-3, j-1} - A_{i, j-1} - A_{i-3,j}$
To convolve with the non-zero submatrix of $H_2$, a matrix in $\mathbb{R}^{3 \times 2}$ you will use $A_{i,j} + A_{i-3, j-2} - A_{i, j-2} - A_{i-3,j}$
Of course you can have other conventions for the indices...


Optimal pulse shape for minimal interference in adjacent frequencies

The root-raised-cosine pulse-shape is widely used because

using the root-raised-cosine as both pulse shape and as matched filter yields the raised-cosine pulse-shape, and
the root-raised-cosine and raised-cosine pulse-shapes are both bandlimited.

The frequency response (Fourier transform) of the root-raised-cosine is equal to the square root of the frequency response of the raised-cosine, so these Fourier transforms are nonzero in exactly the same (bounded) intervals.

Suppose that $h_1(t)$ and $h_2(t)$ are two root-raised-cosine pulses with identical rolloff parameter ($\beta$) whose Fourier transforms are nonzero on non-intersecting intervals. Then matched filtering of of $h_m(t)$ with $h_n(t)$ yields
\begin{equation}
g(\tau) = \int_{-\infty}^{\infty}h_m(t)h_n(\tau - t)dt.
\end{equation}
The Fourier transform of $g$ is
\begin{equation}
G(\omega) ~=~ (\textrm{Fourier transform of $h_m\ast h_n$}) ~=~ H_m(\omega)H_n(\omega).
\end{equation}
Since $H_m$ and $H_n$ are nonzero on non-intersecting intervals, their product is zero for all $\omega$. The inverse Fourier transform of an all-zero Fourier transforms is the zero function, so $g(\tau) = 0$ for all $\tau$. Hence, using a root-raised-cosine matched filter will filter out root-raised-cosine pulses from other frequency bands.

See the answers to this question for more information.


Determine the time length of audio training samples?

It's a very bad idea to stack together consecutive frames in a single vector. By doing so, you are training a classifier to recognize the exact same sound as the one in your training database, pronounced at the exact same speed.

For your application (which appears to be speaker identification), you should consider each individual frame as a training instance (aka "bag-of-frames" models). That is to say, each training sample will give many 12-dimension feature vectors.


what is OQAM (Offset QAM) and why it is required?

Actually, in theory and error performance there is no difference. But in practice, OQPSK  decrases the maximum phase alteriation from 180 to 90 degrees which is better for RF power amplifiers.


The most direct way is to turn your transfer function into frequency response by setting $$s = j\Omega$$ and then $$H(j\Omega) = \frac{j\Omega - 1}{-\Omega^2 + 2j\Omega + 2}$$ So for $\Omega = 0$, you get $$H(0) = -\frac{1}{2}$$ which gives $|H(0)| = \frac{1}{2}$ and is in agreement with your figure. Similarly, for $\Omega = \pm 1$, you get $$H(\pm j1) = \frac{\pm j - 1}{-1 \pm 2j + 2} = \frac{\pm j -1}{1 \pm 2j}$$ You can show that $|H(\pm j1)| = \sqrt{2/5}$, which also agrees with your magnitude response plot.

However, these are only some interesting values of the magnitude response. If you need to draw the magnitude response in greater detail without using some software to do this for you, you have to go for the geometric approach. Here is a good place to start.


Basic Algorithm for Image Classification

Have a look at this report by Marques (PDF). In particular, have a read through section 1.4. It describes several different aspects to the face detection problem. You might want to decide what scenario you're after. Fully unconstrained (uncontrolled environment) face detection is hard to do accurately.

As a first approach, try just doing color-based segmentation.  See this paper by Singh, Chauhan, Vatsa and Singh for some details about how to use color to segment your images for face detection.


What's the difference between "Magnitude Spectrum" and "FFT"?

FFT (fast Fourier transform) is an algorithm that implements the discrete Fourier transform, in a fast manner. Both computes a "full spectrum" from a signal, which is generally complex. The amplitude spectrum is, up to minor corrections (scaling of components), the magnitude of the "full spectrum".


Definition of MIMO (Multi-input Multi-output ),does this output mean antenna or transmitter?

I want to ask the definition of MIMO(Multi-input Multi-output ) If there are one transmitter and one receiver,each of them has N antennas,N>1. So i can say it is a MIMO system.

Any communication system comprising of a transmitter with $N >1$ and receiver comprising of $M>1$ antennas is considered to be a MIMO system. So, this is not only limited to $N=M$ in the way you describe it.

  I think it should be called MISO,because each transmitter just has only one antenna.

No, I think it should be called MU-MISO (MultiUser MISO) because it's what it is; you're transmitting to multiple users (each with a single antenna).

  However,in this paper,it called this system "MIMO" system.

In some applications, MU-MISO and MU-SIMO are seen as MIMO.


Dealing with echoes / dead-spots in digital transmissions

If the channel is static and you are in a deep fade -- then, yes, your only option is to increase the bandwidth or go multicarrier like you suggest.

In fact, modern WiFi systems use that exact technique: OFDM ("Orthogonal frequency division multiplexing" -- sometimes known as multi-tone or multicarrier transmission).

Alternatively, you could go with two receive antennas (microphones in your case) spaced far enough. If you use 2000Hz carrier, the wavelength is v/f = 300[m/s] / 2000 Hz = 0.15 m. So a distance of around 10cm should give you enough diversity. This can also be implemented at the transmitter.


Is it wise to use mathematical morphology and Gabor filter to extract blood vessels in the same algorithm?

In a word, yes.

Mathematical morphology and harmonic analysis (Gabor for instance) are two families of tools, relying on different sets of assumptions, pertaining to image formation models. Harmonic analysis dwells in vector spaces, addressing better linear or superposition  effects. Mathematical morphology, built on lattices, could be more adapted to connected objects, discretized levels, non-linear effects. Both types of contents coexist in images, so mathematical morphology and  harmonic analysis can be complementary tools for complicated tasks.

Some tools reside at the intersection, like morphological wavelets, some links relate linear filters and structuring elements (logarithmically, through the Cramer transform, see fro instance An Explanation for the Logarithmic Connection between Linear and Morphological System Theory). However, it remains wise to use the right set of tools when it is at its best.

For instance, in the most simple signal processing applications, where information is degraded by both  impulse and Gaussian  noise, it is common to first use a non-linear median filter, and a linear low-pass one.

Finally, let me mention a recent thesis related to your topic: RORPO : A morphological framework for curvilinear structure analysis; Application to the filtering and segmentation of blood vessels, where mathematical morphology and harmonic analysis are discussed.


Time-frequency analysis of non-sinusoidal periodic signals

Let us start with the unsupervised methods...

A first approach would be to compute a spectrogram and factorize it with NMF (Non-negative Matrix Factorization). If you are unfamiliar with this technique, it decomposes a spectrogram into a sum of $k$ constant-spectrum sources, each of these having a time-varying amplitude envelope applied to them. This model perfectly suits your problem and there is a good chance that the columns of the decomposition will be spectra for each  waveform/frequency pair (from which you can estimate a fundamental frequency through spectral sum or by turning it back into an autocorrelation function), and the rows will be the activation signals. There are many implementation out there and it's cheap to just throw your problem at it, there's a good chance it'll just work.

Note that the underlying model of NMF is less restrictive than your signal model...

First, because your activation envelopes are either 0 or 1. I haven't dealt with such constraints in the past, but there's probably a way of defining a penalty measure on the terms of the activation matrix $W$ different from 0 and 1, adding that to the optimized criterion and derive a new set of multiplicative update equations. One constraint not mentioned in your comment is that you probably also assume that the sources are not "blinking" in and out rapidly, and stay active/inactive, for a rather long number of consecutive frames. Constraints can be added to NMF to penalize discontinuous activation envelopes. See Virtanen's paper on continuity constraints for that.

The second specificity of your problem is that you might want to make the assumption that all sources have the same waveform, that is to say, that all sources have the same log-frequency spectrum, modulo a translation. To address this specifity, a recommended technique is to compute a constant-Q spectrogram on your input data (which has the right shift-invariance), and perform a shifted-NMF, with one target source. This will recover the spectrum common to all your sources; and the spectral shifts and temporal activation matrix will provide the equivalent of a spectrogram tailored to your source spectrum. See the synthetic signal in the results section of this paper - this is exactly the signal model you describe. This might require some work, but you could probably tweak this method as well to your more rigid signal model in which the activations are either 0 or 1.

Now on to the supervised methods.

If the set of candidate waveforms is known in advance and small and/or the number of frequencies at which they can appear is small (say you want to transcribe music played on a 80s toy keyboard), you could afford matching pursuit. This will correlate your signal with impulse responses consisting of a windowed signal for all searched waveforms/frequencies pair. Computational costs might rapidly make it a bad choice.

Another supervised method is to use the same kind of multiplicative gradient updates as for NMF, but keeping the matrix $W$ containing the spectral templates "locked". That is to say, you build a matrix W with all candidate spectra you want to probe, and perform your multiplicative update to decompose your spectrogram into a product of an activation matrix H and W. Again, update rules tweakable to include the binary activation constraints. See Bertin's work on piano transcription which uses this technique to factorize spectrograms into a sum of synthetic piano note spectrograms (a good start here). This has also been used for drum transcription by Paulus - using drum samples as a basis.


The phenomenon you see is called "window effect" or also "spectral leakage". It occurs because ideally you would like to analyze an infinitely long signal (which would give you ideal peaks in the frequency domain), but actually you analyze a truncated (windowed) signal. In fact you multiply your signal with a window function - in your case a simple rectangular window - which in the frequency domain corresponds to a convolution with the frequency response of the window. This gives you the effect that you see. You can try different window functions, but you will not get rid of the basic effect.

BTW, for the plot you could try using values in dB, i.e. use $20\log (|Y|)$ instead of a linear plot; this will make things look a bit better.


Significance of sequency ordering in Walsh-Hadamard matrices

There isn't much significance. Except for few very particular applications like auto-focusing where you might try to adjust the lens such that the higher sequency components has maximum magnitude.  The FFT of a particular sequency isn't bandwith limited. Therefore sequency has very limited meaning.
You should indeed be looking at the applications like random projections, locality sensitive hashing, associative memory and neural networks:
WHT applications


Turns out I was using same state buffer and instance for all three filters.

Whenever an array finishes filtering, new filter uses the old state array for that different timeline and it results in non-linear outputs at the end/start of a new one.

I created different state arrays and instances for each different array and the problem got fixed. Working smoothly as intended now.


Could a DCT be used for an audio magnitude spectrum rather than DFT?

Yes, DCT can be used to provide a magnitude spectrum with twice the density. I do not quite understand overlap, but I am assuming that since DCT covers less, you thought there would be an overlap. To provide an eligible answer to the question, let me make a quick review for usage of DCT in mainly image processing.

First, we need to make some assumptions. In order to use DCT, you need to have a real signal. This is by definition. While you are saying, DCT has half the bin size comparing to DFT in size N, you are assuming that the signal is low frequency signal. Otherwise, not so much.

For usage of DCT in compression, since DFT of image will be symmetric, it produces redundant information (one side mirror will be enough to reproduce the signal). Therefore, kernel of DCT is used in order to produce denser information comparing to DFT. This is also true for low frequency audio signals, it can be used in the same way. While it makes it denser, coefficients gets bigger, since kernel of DCT covers both sides(real and imaginary parts) of the signal.

My major is image processing, so I tried to map DCT and DFT concepts and explanations in image processing. One difference between image and audio could be sizes, though. In image processing, you know the sizes(row and columns for FFT and other purpose of processing). I guess that you need to divide the vector of audio data somehow in order to further process. Without knowing the data, this could be troublesome(I am not sure).

Here is an image taken from web, but I did not write it down where I took it, could be wikipedia.;

As you can see, transformed image is represented in DCT by magnitude spectrum with no problem. In a more compact and denser way, and look at the magnitude of coefficients. It is bigger than two times of DFT. DFT is symmetric, you could just divide it into two. One part is redundant. And one more thing, DCT can store the information is not just half of DFT but nearly quarter of DFT. That is generally the case of DCT overcoming to DFT in images.


How do I know quantitatively if the correlation of two time series is significant?

I believe you are looking for the normalized cross-correlation coefficient (source), which in your case would be computed as follows. First, you de-mean your time series $x(t)$ and $y(t)$ such that

$f(t) = x(t) - \bar{x}$, and
$g(t) = y(t) - \bar{y}$

You can treat your time series as vectors $\mathbf{f}$ and $\mathbf{g}$. Assuming the vectors are of equal length, the formula for the normalized correlation coefficient is:

$R_{\text norm}  = \frac{\mathbf{f}^T \mathbf{g}}{\| \mathbf{f} \| \| \mathbf{g} \|}$

This is bound to have a maximum value of $1$ (meaning the time series are very correlated), as you are essentially multiplying together two unit vectors.

However, note that if you have different length input vectors, then you should not scale the output.


The relationship between downsampling and frequency resolution

Your understanding is correct. I don't have the full text of the paper, but it sounds like they aren't being very precise with their description. As you pointed out, downsampling doesn't improve frequency resolution in and of itself, it merely reduces the amount of bandwidth that can be represented unambiguously.

However, I can take a guess at what they were referring to. Assuming they are using a DFT for spectral analysis (the common choice), note that the fact that for a given DFT size $N$, the output bin spacing is a function of the sample rate:

$$
f_{res} = \frac{f_s}{N}
$$

So, if a smaller frequency resolution in your spectral calculation is desired, decreasing the sample rate of the data can help you do that. Of course, you could also increase the DFT size $N$ as well; there are some applications where that's not an option, though (e.g. if you're using hardware that is limited to certain FFT sizes). Also, this method can have advantages when the signal of interest is in a small band of frequencies relative to the original sample rate; in this case, downsampling followed by a smaller DFT can be more computationally advantageous than doing one large DFT (this is known as the zoom FFT technique).

Note that there isn't a free lunch here: to actually realize an increase in resolution, you need $N$ samples at the lower sample rate, which requires a longer observation time. This is the same total observation time that is required by just collecting more data at the higher sample rate and then doing a larger DFT.


linear phase and generalised linear phase filters

Yes. The time delay of real-coefficient linear-phase N-point FIR filter is (N-1)/2 samples. The time delay of complex-coefficient generalized-linear-phase N-point FIR filter is also (N-1)/2 samples.

You can prove this to yourself.  Design a narrowband linear-phase lowpass FIR filter, and plot its group delay. Then multiply that filter's coefficients by a complex exponential whose frequency is one fourth the sample rate, creating a complex-coefficient filter. Now plot the group delay of the complex filter. Compare the two group delay plots.


Question regarding 16 QAM digital modulation

First of all, I recommend getting familiar with Matlab's doc command. doc qammod has the answers to your second question: qammod requires a vector of integers between 0 and M-1.

Regarding your first question: the reason to reshape your bits into groups of four is that 16-QAM transmits $\log_2(16)=4$ bits per symbol. So, your bi2de command requires binary numbers of four bits to get decimal integers betwen 0 and 15.

If you want to save lines of code, what you could do is generate integers between 0 and 15 directly, with dataSymbolsIn = randi([0 15], n/4, 1). However, doing this has one disadvantage: calculating the bit error rate at the receiver is more complicated if all you have are the integers. Let's say that you transmit a 15 and receive a 13 (obtained with qamdemod), how many bit errors there were? This is why it's common practice to start with actual bits in the transmitter and compare them one-to-one with the received bits.


calculate fft convolution for an LTI system given its input (in time domain) and its poles and zeros

In general the filter you want to implement is an infinite impulse response (IIR) filter, unless all poles are at the origin of the complex plane (assuming causality), in which case it is a finite impulse response (FIR) filter.
In the general (IIR) case, your suggested method will not result in an exact implementation of the filter. There are two reasons for this:

by using an FFT of the filter's frequency response you actually approximate the infinitely long impulse response by an impulse response of finite length. That impulse response is simply given by the inverse FFT of the sampled frequency response.

multiplication of the FFTs implements circular (cyclic) convolution, which is different from linear convolution.

Both errors can be made small. For the first error (FIR approximation of IIR filter), you just need to choose an FFT length that captures a large percentage of the energy of the impulse response. This basically means that you choose an FIR filter of sufficient length to approximate the given IIR filter. For minimizing the second error (circular convolution instead of linear convolution), you need to zero-pad the input sequence and the FIR approximation of the IIR filter.
The question remains why one would want to use an implementation like this. One disadvantage of the proposed solution is that you have to wait for the complete input signal before you can start processing, i.e., you introduce a substantial delay. This problem could be tackled by block processing, like overlap-save or overlap-add in the case of standard fast convolution FIR filtering. The other disadvantage is the increase in memory requirements and computational load. Most practical IIR filters have relatively low orders (lower than $20$), but an FIR filter that provides a reasonable approximation will usually have many hundreds of coefficients, or even more.
There exists an exact method for block processing of IIR filters, in which FFTs can be used to solve certain matrix-vector multiplications. This method is explained in detail in this very accessible book chapter by Selesnick and Burrus: Fast Convolution and Filtering (Section 8.3.1).


Is there any tool to anlyze audio specturm of a song and give the average visually ?

Audacity will let you look at a song's spectrum/spectrogram, and it works with most audio formats.

It's also free!


There are errors in the presented answer. Here's the same but corrected and rewritten a bit:

 1.01 * 1.11
--------------
 carry: 1 1
        1.1111  (shift right 1.11 with sign extension)
        0.000
      + 0.01    (2's complement of 1.11)
--------------
      = 0.0011

Here, 2's complement means negation, which is done by flipping each bit and adding one least significant bit valued 1. In multiplication of signed binary numbers, if the first multiplicand (the one that is not shifted) is negative, the correct result is obtained by doing subtraction in place of the last addition by the second multiplicand, and that's the + 0.01 line. It's a binary trick that is equivalent to actually dragging along all the sign extensions, which leads to large values in the carry and more steps than necessary:

 1.01 * 1.11
--------------
 carry: 2 1
        1.1111  (shift right 1.11 with sign extension)
        0.000
        1.11
        1.1     (sign extension of the multiplicand 1.01 is ...111.01)
      + 1.
--------------
      = 0.0011

There is no need to have two binary decimals before . in the result.

The subtraction trick can be tested with different values and it passes all the tests:

              carry:             1               1
second multiplicand: 000 001 010 011 100 101 110 111
                     00  01  10  11  00  01  10  11
                   + 0   1   0   1   0   1   0   1
----------------------------------------------------
                   = 000 111 110 101 100 011 010 001
     2's complement: 000 111 110 101 100 011 010 001


Is Wavelet Transform better than FT if I'm only interested in knowing the frequencies and not interested of when did they happen?

Since they seem more fashionable than Fourier transforms, some persons tend to overuse wavelets in contexts where they do not belong, or at least don't do a much better job.

For "knowing the frequencies", Fourier transforms are the correct tools to use first: it is so easy to badly select a wavelet transformation (including its shape and its sampling) that will give you outcomes difficult to interpret, and even erroneous information. To be more specific, among fast orthogonal transforms (DWT and FFT), due to the restrictions on dyadic wavelets, it is unlikely to have better results with wavelets than Fourier (because of shift-variance, aliasing and badly-shaped FIR wavelets).

  does WT do any better job than FT in detecting the frequencies

It really depends on what you expect in detection, and the nature of your signals, and it often takes more than a transform. However, I would recommend to use wavelets or time-frequency tools if Fourier reveals not sufficient, or in parallel to address statonnarity, etc., but  not instead.


You complex multiplication code inputs 2 samples for every output sample.  That will act as a frequency doubler.


Direct Inverse Filtering an image, spatially convolved with a blurring function without any noise, doees not give my original image back?

I think I found out the answer to my own question.

Since I'm using an averaging filter mask (some constant multiplied by a matrix of ones), the Fourier Transform of such a filter is an impulse at zero frequency.

When I pad the mask filter in step 1 to make the dimensions of the mask equal to the image, I am 'truncating' the filter, ie. windowing the mask, which is akin to multiplying the Fourier Transform of the mask with a sinc function.

As a result, the Frequency Domain filter I use in step 5, ie. H, is not an impulse but a sinc function. That's the reason why my results are not ideal.

On the other hand, if I were to use an impulse in step 5, only the zero-frequency value would remain, and my result would be a constant image.

In short, it's a lose-lose situation.


Rendering a spectrogram of audio data - resulting magnitudes are all too close together

A periodic looking output from an FFT could mean a bug such that the FFT did nothing (output == input, unchanged), or that the input to the FFT only had a single large magnitude or non-zero data point (input == garbage).


I would do it like this:

Start with a 2 MHz signal containing up to 10 FM stations. In practice it will contain maybe three or so.
Identify the carrier frequencies of each station; call them f0 to f9. These are the frequencies in your downconverted 2 MHz signal, not the original carriers.
Design 10 200kHz-wide bandpass filters centered around each f0:f9. Filter the 2-MHz signal with each of these filters to separate each individual FM carrier.
Downconvert each FM carrier back to baseband. To do this, multiply each carrier by the signal exp(-j*2*pi*fi*t), where fi is the carrier frequency and t is time. This will bring each FM signal to baseband.
Demodulate each signal using a regular FM demodulator.
Send each resulting audio signal to its destination.

An alternative procedure would swap steps 3 and 4:

Frequency-shift the entire 2 MHz signal so that the desired carrier ends up at 0 Hz.
Apply a low-pass filter to remove all undesired frequencies.

The advantage of this method is that you need to design only one low-pass filter instead of 10 band-pass filters. However, you need to be careful that you don't introduce aliases that can affect the frequencies you want to recover. The reason is that a 2-MHz signal sampled at 4-MHz doesn't have any "headroom" at the ends; any time shifting will introduce aliases. Just make sure these don't end up in the passband of your low-pass filter.

So, you need to: design bandpass filters (or one low-pass filter), generate complex exponentials, and demodulate the FM signals. You can do this Matlab, Liquid DSP, GnuRadio, or any number of other numerical/radio programs and languages. I think your largest constraint will be computer power, unless you only have a couple of FM stations.


Amplitude problem by puting 5 FIR in a row

Look at the power spectral density of the signal prior to application of filter b5.

If you filtered white noise with two filters, one with a bandwidth of 11 MHz and the other with a bandwidth of 60 MHz, and both with the same passband gain, then indeed the difference between the power of the two output signals should be about 10*log10(60/11) or 7.3 dB.  But if, for example, your signal was already bandlimited to the 11 MHz bandwidth of filter 1, then there would be little difference between the output powers of the two filters.  A signal with a more complicated spectra (e.g., signal is within 11 MHz, but broadband noise exists) will be somewhere in between.

Also, if you are trying to compute the power of the signal in dB, you would not have a sqrt in your 10.*log10(sqrt(I.*I + Q.*Q)) equation.


The convolution in the time domain is translated to simple multiplication in the frequency domain. This means that since the values of the transform of x(t) are at the integer locations, the result of the multiplication will be the transform of x(t), as you defined the integer values for h(t) transform to be 1 at integer bins.


UPDATE: I suspected Spectral Leakage but that is not the case here. The reason for this is the OP is choosing to plot just the real component of the FFT, with a linear magnitude scale so the actual effect of the spectral leakage cannot be seen. (We actually get spectral leakage when we use endpoint = true!). If you change the plot from fplt.real to np.abs(fplt) you will see the true magnitude plot that must include the real and imaginary components!).

Here is the result with endpoint = true and plotting the dB of the magnitude.
plt.plot(tplt, 20*np.log10(np.abs(fplt)), 'b-'). The spectral leakage in this case is quite clear!

Here is the result with endpoint = false using
plt.plot(tplt, 20*np.log10(np.abs(fplt)), 'b-')

Here we see noise that is nearly 300 dB down (floating point precision).


How is "Joint class histogram" algorithm deploy using Matlab?

I get the answer from another source and I deploy like this

% load floating and fix image

%% image segmentation using k-mean clustring.
idx_fix = kmeans(double(fix(:)),4);
idx_float = kmeans(double(float(:)),4);

%% display image segmentation
mat1 = reshape(idx_fix,[],500);%besed on your image size
mat2  = reshape(idx_float,[],500);
imshow(mat1/max(mat1(:))); title('Fix image class');
figure;
imshow(mat2/max(mat2(:))); title('Float image class');

%% Joint class histogram
joint = zeros(5,5);%besed on number of class
for i=1:length(idx_fix)
   c1 = idx_fix(i) + 1;%offset matrix by 1 because index start by 1
   c2 = idx_float(i) + 1;
   joint(c1, c2) =  joint(c1,c2) + 1;
end

If the image aligned correctly the joint class histogram result as a diagonal line.

Thanks, Walter Roberson


Instead of having a hard edge between the image data of interest and the background, you could use a two-dimensional tapered window function, as is often done in spectral analysis. You might start by trying a Gaussian window, which for a two-dimensional case would look something like:

$$
w[x,y] = e^{-\frac{\left(x-\frac{N_x-1}{2}\right)^2}{2\left(\sigma_x \frac{N_x-1}{2}\right)^2}} e^{-\frac{\left(y-\frac{N_y-1}{2}\right)^2}{2\left(\sigma_y \frac{N_y-1}{2}\right)^2}}
$$

$N_x$ and $N_y$ are the dimensions of the desired transform in the $x$ and $y$ directions, respectively, and $\sigma_x$ and $\sigma_y$ are parameters that allow you to control the shape of the window; for small $\sigma$ values, most of the energy in the window function will be concentrated toward the center, with that effect decreasing as you increase $\sigma$.


ARMA models are useful when you need to model a Signal plus Noise situation where the signal is an AR process and the noise models sensor noise. The overall model is an ARMA model.

See HL Van Trees, Detection, Estimation, and Modulation Theory, vol 4 Array Processing. He gives an example of a Spatial AR process sensed by noisy sensors. The overall model is ARMA

ARMA model are also useful in situations where you have strong nulls n your spectrum like multi path and you have a noise like signal


If you have complex conjugate poles (i.e., $\text{Im}\{p_i\}\neq 0$) you have oscillations in the step response. Only real-valued poles will give you a monotonic step response as shown in the figure.

Take a look at this related answer.


I need advice about how to make an audio frequency analyzer

This really depends on how fancy you want to get. A "good" analyzer will typically do the following things.

Split the input into frames. The frames typically overlap and are windowed. Good choices are an overlap of 50% and a Hanning window.
Do an FFT
Select the center frequencies of bands. For audio good choices are octaves or third octaves (see http://www.engineeringtoolbox.com/octave-bands-frequency-limits-d_1602.html)
For each integrate the FFT energy around the center frequency. There are different choices for the integration windows. Good choices are trapezoids or butterworth windows. A really good window choice will make sure that all the energy is accounted for, i.e. that the sum of the band energies equals the sum of all FFT bins.

If the number of frequency bands is relatively small, it is more efficient and easier to simply use a bank of parallel bandpass filters and an RMS detector. This also allows to dial in the temporal behavior of the display nicely.


How do you show the Discrete Time Fourier Transform of $x[n]=\cos(2\pi f_0n)$ is $ \frac{1}{2} \delta(f+f_0) + \frac{1}{2}\delta(f-f_0)$?

The sequence $ ~x[n] = \cos(\omega_0 n) ~$ , $-\infty < n < \infty$, is neither absolutely nor square summable, therefore its DTFT formally does not exist; i.e., the DTFT sum does not converge to a finite number, but diverges to infinity.
Because of the extreme importance of that signal in the context of both the theory and the practice of signal processing, we would like to express its DTFT in an adequate and consistent manner under all operations of signal processing. One such possible representation makes use of the Dirac impulse functions with the following observation:
Let $x[n] = 1$, for all $n$, be a discrete-time (unstable) signal whose DTFT  $X(\omega)$ given by the sum: $$ X(\omega) = \sum_{n=-\infty}^{\infty} 1 e^{-j \omega n} \tag{1} $$
This sum, resembles a CTFS (continous time Fourier series) synthesis relation, remembering the fact that the DTFT $X(\omega)$ is a periodic function of continuous -frequency $\omega$. Then the above sum in Eq.1 is re-interpreted as a CTFS synthesis as:
$$ X(\omega) = \sum_{n=-\infty}^{\infty} a_n e^{-j \frac{2\pi}{W} n \omega } \tag{2} $$ where  $\frac{2\pi}{W}$ is the fundamental frequency with the fundamental period being $W = 2\pi$, and the CTFS coefficients are all $a_n = 1$. From CTFS pairs, we find the signal $X(\omega)$ whose CTFS coefficents are all $1$, as the Dirac impulse train :
$$ X(\omega) = 2\pi \sum_{n=-\infty}^{\infty} \delta(\omega - 2\pi n) \tag{3}$$
and using the CTFS analysis equation, we have :
$$
\begin{align}
a_n &= \frac{1}{2\pi} \int_{-\pi}^{\pi} \left( 2\pi \sum_{n=-\infty}^{\infty} \delta(\omega - 2\pi n) \right) e^{-j \frac{2 \pi}{2\pi} \omega n } d\omega \\\\  &= \frac{1}{2\pi} \int_{-\pi}^{\pi} 2\pi \delta(\omega) e^{-j \omega n } d\omega = 1 \tag{4}
\end{align}
$$
The integral in Eq.4 is identical to an inverse DTFT of the impulse function $2 \pi \delta(\omega)$, hence we deduce that the DTFT $X(\omega)$ of $x[n]=1$ is
$$ x[n] = 1 \implies X(\omega) = 2\pi \delta(\omega) \tag{5}$$
Having found the DTFT of the sequence $x[n]=1$, we apply an estabslihed property of the DTFT :
$$ x[n] \leftrightarrow X(\omega) \implies x[n]e^{j\omega_0 n} \leftrightarrow X(\omega - \omega_0) \tag{6}$$
$$ 1 \leftrightarrow 2 \pi \delta(\omega) \implies 1 \cdot e^{j\omega_0 n} \leftrightarrow 2\pi \delta(\omega - \omega_0) \tag{7} \\\\$$
to the seqeunce $x[n] = \cos(\omega_0 n)$ after decomposing it into complex exponentials:
$$
\begin{align}
x[n] &= \cos(\omega_0 n) = 0.5 \left( e^{j \omega_0 n} + e^{-j \omega_0 n} \right)  \\ \\
X(\omega) &= \pi \delta( \omega - \omega_0) + \pi  \delta(\omega + \omega_0)
\tag{8} \\ \\
\end{align}
$$
Eq.8 gives us the requested DTFT of the sequence $x[n] = \cos(\omega_0 n)$.


How to calculate resolution of DFT with Hamming/Hann window?

The frequency resolution of a DFT is indeed significantly affected by the choice of window. Two key parameters to consider when selecting a window are the frequency resolution and dynamic range. The frequency resolution is the ability to discern two closely spaced frequencies at similar power levels (therefore a metric of what that spacing is would be the frequency resolution), while dynamic range is the ability to discern two signals at different frequencies and significantly different power levels; specifically the ability to see a low level signal in the presence of a much stronger signal (therefore a metric of what this power difference is would be the dynamic range). Typically a higher dynamic range is achieved by using appropriate windowing at the expense of frequency resolution. The highest frequency resolution that can be achieved is with a rectangular window (and the lowest loss; any other window will also introduce a processing loss).

The Equivalent Noise Bandwidth (ENBW) is a commonly used metric for frequency resolution. The ENBW is the bandwidth (often given in FFT bins) of a brickwall filter that would result in the same noise power as the DFT "filter" for white noise (the DFT can be described as a bank of filters). For the rectangular window, the ENBW is one bin, but the poorest dynamic range (sidelobes with peaks rolling off at a $1/f$ rate in magnitude.  Any other windowing used would have a higher ENBW.

For a given length N window function $w[n]$, the ENBW in bins is computed as follows:

$$ENBW = N \frac{\sum \left(w[n]\right)^2 }{\left(\sum w[n] \right)^2}$$

See this classic paper by fred harris tabulating the ENBW as well as other paramters of many commonly used windows. fred harris windowing

Below is a slide I have further depicting the ENBW for a rectangular window:

Below compares a N = 30 rectangular window (denoted the Dirichlet kernel) to a Blackman window of the same length. The decrease in level shown of -7.83 dB is the coherent gain for this particular window used. The coherent gain is simply the average of the window weights. The result below is is consistent with the coherent gain of 0.42 as listed for the Blackman window tabulated in the fred harris reference given; $20Log_{10}(0.42) = -7.54$ dB.  What is very clear from this picture is the significant increase in both dynamic range and frequency resolution.

The reason for not being an exact match to the tabulated results is the window length. As the window length approaches infinity the physical result approaches the tabulated, as shown in the figure below for coherent gain. For this reason caution must be applied when using the tabulated results for any smaller length windows.


3.

  The polar image has the black areas, which are parts that can't be
  mapped to the original image because it would be outside of its
  boundaries. However, these areas have nearly disappeared in the log
  polar image. Why is that?

Look at the second (polar) and third (log-polar) images.  What's different about them?

In the 3rd image, stuff that was on the left is now in the center.  Stuff that was in the center is now on the right.  The "log" in log-polar means that you're stretching the image in a logarithmic way.  So the black out-of-image scallops are still there, they've just been squished very thin against the right side of the image, while details from the left side of the image have been stretched out to take their place.


effect of high pass filtering and low pass filtering on background color

It's very simple.  White = 255, Black = 0.  After high pass filtering, the pixels that don't change relative to their neighbors (low spatial frequency) get filtered out--meaning the pixel value is reduced.  Thus, any solid color areas will be replaced by blackness, and only points of pixel value transition will be white.


Phase tracking in communication receivers

The main purpose of a phase synchronization system, whether Costas loop or a PLL, is to track the phase of an input signal. The bandwidth of a PLL determines whether the signal falls within the range it can track or suppress. If the signal can pass through and the bandwidth is narrow, less noise enters the system increasing the loop SNR, while a wider bandwidth allows more noise deteriorating the tracking performance. But do remember that the narrower the bandwidth, the longer the acquisition time due to the smaller variation injected by the PLL.


What is the relationship between image spatial domain and FFT powerspectrum?

The Contrast Transfer Function (CTF) is essentially equivalent to the frequency response of your system.

Therefore, your question is about how to read the two dimensional Discrete Fourier Transform and how to interpret its content.

To relate what is mentioned in that answer to your specific question, consider the following CTF, obtained from Wikipedia (Axis added by me):

This is a two dimensional spectrum. If you are familiar with a one dimensional spectrum, it is as if you took that one and rotated it about the centre of the image.

Therefore, low (spatial) frequencies are towards the centre of the image and high (spatial) frequencies are towards the edges of the image.

In fact, if we take a cross section of that image (along the direction of the blue line), we see something like this:

And if we shift the origin of the spectrum to its conventional position (for a 1D spectrum), we get:

So, this now describes how well each spatial frequency goes through the optical system assuming, it is symmetrical.

To relate this to physical frequency, we have to do a bit more work. Now, this is a simple grayscale image, which means that its maximum value is 255. So to normalise the "height" of this curve, we need to simply divide by 255.

The other thing we have to do is to normalise the frequency. This image is 512x512, the DC (or, constant value) is at 0, the Fs is at 512 and the Nyquist frequency is at 256. The formula is a plain simple division of $f = \frac{k}{XDim} \cdot Fs$ where $k$ is the frequency "bin" (of the discrete spectrum) and the fractional part of this we can simply represent in the graph.

Taking all of this into account gives us:

So, now, whatever your Fs turns out to be, you can put it on the diagram and work backwards to find how much each spatial frequency gets diminished as it passes through this system.

Spatial resolution is measured in "lines" per unit of length. In your case, your pixel size is 1.7Å. The "highest" frequency (in any direction) that can theoretically go through a system that is digitised at that level is one black line of width 1.7Å followed immediately by another white line of width 1.7Å and that is one cycle of period 3.4Å.

Therefore, your Nyquist frequency is 1 line per 3.4Å or $\approx 2941176471$ lines per mm.

Hope this helps.

EDIT:

All the plots in this post were created using Octave and the following code:

% Load the image processing package
pkg load image
% Load the image that was downloaded from Wikipedia
I = imread("~/Downloads/CTF.jpg");
% The image is 512x512, let's take the cross section of the image to see a slice of the spectrum
plot(I(256,:));
% Now shift the origin using fftshift
plot(fftshift(I(256,:)));
% Add normalised X and Y axis, labels, grid, etc
plot((0:511)./512,double(fftshift(I(256,:)))/255.0);xlabel("Normalised Spatial Frequency");ylabel("Normalised Transfer Coefficient");grid on;


Constant Modulus Algorithm and the Gradient Operation

As mentioned in the comments, the symbol $^H$ denotes the conjugate transpose of a matrix or vector, which means that the vector/matrix is transposed and that all of its elements are conjugated. The bar over $y_k$ means complex conjugate (note that $y_k$ is a scalar, not a vector or matrix).

When computing the gradient with respect to a complex variable in order to maximize or minimize a function, there's a trick which is explained in more detail in this answer. You basically take the derivative with respect to the conjugated variable and regard the non-conjugated variable as a constant. So when computing the gradient of the expression $\mathbf{w}^H\mathbf{x}\mathbf{x}^H\mathbf{w}$, you formally take the derivative w.r.t. $\mathbf{w}^H$, regarding $\mathbf{x}\mathbf{x}^H\mathbf{w}$ as constant, which yields the given result.

As for your third question, you should realize what you actually want to do. You want to minimize the error function, so you do not want to use a steepest ascent but a steepest descent algorithm, moving in the direction of the negative gradient.

And the Matlab code does in fact implement the equations you stated (replacing the modulus by the variable R2, which is equal to $1$ in your equations). Just use the error $e_k=|y_k|^2-1$ and the expression $y_k=\mathbf{w}^{(k)H}\mathbf{x}_k$ in your update equation, and you'll see it.


Periodicity of the discrete-time Fourier Transform

The DTFT is always $2\pi$-periodic. However, it can also have a smaller period, namely a fraction of $2\pi$. Take any sequence $x[n]$ for which the DTFT exists and insert $L-1$ zeros between the samples. The DTFT of the new sequency $\hat{x}[n]$ can then be written as

$$\hat{X}(e^{j\omega})=\sum_{n=-\infty}^{\infty}\hat{x}[n]e^{-jn\omega}=\sum_{n=-\infty}^{\infty}\hat{x}[nL]e^{-jnL\omega}\tag{1}$$

$\hat{X}(e^{j\omega})$ as given by $(1)$ clearly has a period of $2\pi/L$.

Concerning your second question, the reason why a continuous-time Fourier transform $X(j\omega)$ cannot be written as $X(e^{j\omega})$ is that

$$X(j\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\neq \int_{-\infty}^{\infty}x(t)(e^{j\omega})^{-t}dt\tag{2}$$

simply because generally

$$e^{-j\omega t}\neq (e^{j\omega})^{-t}\tag{3}$$

unless $t$ is an integer (as is the case with the DTFT).

Note that if $(3)$ were not true,  i.e., if $e^{-j\omega t}= (e^{j\omega})^{-t}$ were true, we could easily show that $e^{-j\omega t}=1$ for all $t$. Just write $\omega=2\pi/T$ and you get

$$e^{-j\omega t}=e^{-j2\pi t/T}=\left(e^{-j2\pi}\right)^{t/T}=1^{t/T}=1$$

which is of course absurd.

Proof of Eq. $(3)$:

Let $z$ be a complex number with $|z|=1$ (the magnitude is irrelevant here):

$$z=e^{j\theta}\tag{4}$$

Let $a$ and $b$ be real numbers. Then

$$z^{ab}=e^{j\theta a b}\tag{5}$$

And with $z^a=u$

$$\left(z^a\right)^b=u^b=e^{j\arg\{u\}b}\tag{6}$$

With $\arg\{u\}=\text{pv}\{\theta a\}\in (-\pi,\pi]$, where $\text{pv}$ denotes the principal value, it is clear that generally

$$z^{ab}=e^{j\theta a b}\neq \left(z^a\right)^b=e^{j\,\text{pv}\{\theta a\}b}\tag{7}$$

There are two cases where $(5)$ and $(6)$ are equal:

if $\text{pv}\{\theta a\}=\theta a$, which is the case if $\theta a\in (-\pi,\pi]$.
if $b$ is integer, since

$$\text{pv}\{\theta a\}=\theta a+2\pi k\tag{8}$$

with some appropriately chosen integer $k$ (such that the result is in the interval $(-\pi,\pi]$), then

$$\left(z^a\right)^b=e^{j\,\text{pv}\{\theta a\}b}=e^{j(\theta a+2\pi k)b}=e^{j\theta ab}=z^{ab},\qquad b\in\mathbb{Z}\tag{9}$$

[Note that $(9)$ would also hold for rational $b$ as long as $kb$ is integer.]


If you have a $T$-periodic function $g(t)$ with Fourier coefficients $c_n$

$$g(t)=\sum_{n=-\infty}^{\infty}c_ne^{in\omega_0t},\quad\omega_0=\frac{2\pi}{T}$$

and you create a new function $f(t)$ with Fourier coefficients $d_n$ defined by

$$d_n=c_{n/2},\quad n\textrm{ even}\\
d_n=0,\quad n\textrm{ odd}$$

i.e. the spectrum is expanded by a factor two by inserting zeros between the original Fourier coefficients. Then $f(t)$ is given by

$$f(t)=\sum_{n=-\infty}^{\infty}d_ne^{in\omega_0t}=
\sum_{n\textrm{ even}}c_{n/2}e^{in\omega_0t}=
\sum_{n=-\infty}^{\infty}c_ne^{i2n\omega_0t}=g(2t)$$

So you get a compressed time domain signal with period $T/2$, as directly shown in the first equation of the problem formulation in your question.


Is there any way of retrieving phase information with only access to real part of successive FFT frames?

webaudio isn't returning the real part of the complex FFT transform, but just the log magnitude of the FFT result.  From just this result alone, of the fluctuations between the results of multiple frames, you can't retrieve the phase information.

However, if you can get at the input data and decompose this input vector into a symmetric vector plus an anti-symmetric vector, the magnitude of the FFT of the symmetric vector will be equivalent to scaled real components of the FFT of the original input vector, and the magnitude of the FFT of the anti-symmetric vector will be equal scaled imaginary components of the FFT of the original input vector.  Thus, the decomposition plus 2 magnitude-only FFTs will allow you to determine both the magnitude and phase of the FFT components.

v_even = v + v_reversed
v_odd  = v - v_reversed


How to smooth (or interpolate) phase of FFT and reduce data points

To downsample an FFT result for both magnitude and phase, it may work best to do an FFTShift before the FFT, then downsample the real component vector and then the imaginary component vector separately, before taking the atan2() of them to estimate downsampled phase.


Determine sequence's value from equally null circular convolution DFT

Your professor is wrong, the solution is not unique. There are two (real-valued) solutions:

x1 = [2, 2, 0, 0]; y1 = [1, -1, 1, -1]

with DFTs

X1 = [4, 2-2i, 0, 2+2i]; Y1 = [0, 0, 4, 0]

(the product of which is obviously zero)

and

x2 = [2, -2, 0, 0]; y2 = [1, 1, 1, 1]

with DFTs

X2 = [0, 2+2i, 4, 2-2i]; Y2 = [4, 0, 0, 0]

(the product of which is also zero).

This is how you obtain these two solutions: first of all, note that the DFT of the circular convolution of the two sequences is equal to the product of the DFTs of the two sequences, where the individual DFTs are of length $4$. The product is zero for all indices $0\le k<4$ if for each index at least one of the two DFTs is zero. Since I assume that we're talking about real-valued sequences, we only need to consider the indices up to and including Nyquist (i.e. index $k=2$), because the last index $k=3$ is redundant due to the symmetry of the DFTs.

Since none of the two sequences can be identical to zero ($x[0]$ and $y[0]$ are both given and non-zero), only specific DFT bins can be zero, not all of them. For the length $4$ DFT of the zero-padded length $2$ sequence $x[n]$ it is easy to see that its DFT is given by

$$\begin{align}X[0]&=x[0]+x[1]\\
X[1]&=x[0]-jx[1]\\
X[2]&=x[0]-x[1]\\
X[3]&=x[0]+jx[1]\;\;(=X^*[1]\;)\end{align}\tag{1}$$

from which it is clear that with $x[0]\neq 0$, only $X[0]$ or $X[2]$ (but not both) can be zero. $X[1]$ can never be zero (for real-valued $x[1]$), from which it follows that $Y[1]$ must be zero. Since $Y[k]$ cannot be zero for all $k$, $X[0]$ or $X[2]$ must be zero. This leaves us with two possible choices for $x[1]$:

$$x[1]=x[0]\rightarrow X[2]=0\tag{a}$$

or

$$x[1]=-x[0]\rightarrow X[0]=0\tag{b}$$

In the first case, we need a DFT $Y[k]$ which is zero for all $k$ except for $k=2$, in the other case $Y[k]$ must be zero for all $k$ except for $k=0$. Note that $k=0$ corresponds to DC, and $k=2$ corresponds to Nyquist. The DFT of a sequence is zero everywhere except at DC if the sequence is constant, and it is zero everywhere except at Nyquist if the sequence is alternating. So for choice $(a)$ we need to choose $y[n]$ to alternate, and for choice $(b)$, $y[n]$ needs to be constant. These are the only two real-valued options.

EDIT: If you allow the sequences $x[n]$ and $y[n]$ to be complex-valued, there are two additional solutions. From $(1)$ it can be seen that for $x[0]\neq 0$, $X[k]$ can only be made zero for one index $k$. This means that $Y[k]$ must be zero for all indices except for that one index where $X[k]$ is zero. Therefore, $Y[k]$ must be a spectral line, i.e. its corresponding sequence $y[n]$ must be a complex exponential on the DFT ($N=4$) grid:

$$y[n]=e^{j\pi nl/2},\quad l=0,1,2,3\tag{2}$$

For $y[n]$ as given by $(2)$, $Y[k]$ is zero for $k=l$. The solutions $l=0$ (constant sequence), and $l=2$ alternating sequence are real-valued and were given above. The two other solutions are

y3 = [1, i, -1, -i]

($l=1$) and

y4 = [1, -i, -1, i]

($l=3$). The DFT of the corresponding sequences $x[n]$ must have a zero at $k=l$. From $(1)$ we get

x3 = [2, -2i, 0, 0]
and
x4 = [2, 2i, 0, 0]

The DFTs are

X3 = [2 - 2i, 0, 2 + 2i, 4]; Y3 = [0 4 0 0]

and

X4 = [2 + 2i, 4, 2 - 2i, 0]; Y4 = [0 0 0 4]

Both products are clearly zero.

In sum, there are four solutions to the problem, two of which are real-valued.


If my guess is right, your code show you trying remove the vocal from  second audio, the first audio must be the original vocal and the second is some mix (vocal + music)...

So if my guess is true start trying a simple approach to test, here my simple test.

I get just 4096 points of some audio voice and mix with a sine 4096 points in 200hz sampled in 1600hz, all mono.

Plot 4096 voice points:

Now plot a sine in 200 Hz:

Now the Mix (sine + voice):

Now to get the original sine without the voice you need do:

real((ifft(fft(mix) - fft(voice_Signal))))

The final plot after the step above:

course is a very simple test, and in this case we're telling you to have the original clean audio voice before being mixed ...


Spatial correlation is relevant in sensor arrays and MIMO wireless communications. In the case of MIMO communications where multiple spatially separate antennas receive different observations of a multipath signal, spatial correlation can be used to evaluate the 'antennae diversity'. Higher the diversity, better the mitigation of signal fading effects introduced by the wireless medium. As a result the link quality is improved.
http://en.wikipedia.org/wiki/Antenna_diversity


Plotting the frequency response of a FIR filter

If you want to do it manually, it is not upsampled by adding zeroes between coefficients. It is usually done by adding extra zeroes after the real coefficients. When padded to some convenient power of two like 512 coefficients, FFT can be used. The resulting amplitudes need to be scaled due to the padding. Matlab and Octave have freqz function to plot this directly for you.


Is DFE equalizer always blind or can we use reference signal?

I can answer for adaptive filter containing feedback path alone as I have implemented it recently (both Blind DFE and with Reference Symbols), but it should work for DFE containing feed-forward path as well. In Blind DFE while computing the error, $e[n] = x[n] - \hat{x}[n]$, where $x[n]$ is the equalizer output, $\hat{x}[n]$ is the corresponding Hard Decision. But if you have reference symbol $r[n]$, you need not compare with hard decision output and $e[n] = x[n] - r[n]$.

If the reference symbols are placed in the beginning of data, then the filter adapts more reliably (delay and tap values) in the beginning itself, you can then use your Blind DFE for the remaining data part.

The example I have used comes from 802.11ad Single Carrier data packet where for each frame containing 512 symbols, the first 64 are known reference symbols and the rest 448 are random data symbols (BPSK/QPSK/16QAM). So for the first 64 symbols, I skipped the hard decision slicer and directly compared adaptive filter output with reference symbols. For the remaining 448 symbols, I used Blind DFE.


Relationship between signal rate, data rate, minimum bandwidth

Here's an attempt at an answer, by comparing your equations with the more commonly used ones.

We can write the "pulse rate" or the "signaling rate" $R_p$, measured in baud, or equivalently in pulses per second or symbols per second, as: $$R_p = \frac{1}{k} R_b,$$ where $R_b$ is the bit rate, and $k$ is the number of bits carried by a pulse.

To elaborate: a digital communications signal can be written as a sequence of pulses of different amplitudes, each pulse time-shifted an amount $T_p = 1/R_p$, and with the information carried in the pulses' amplitude: $$s(t) = \sum_j a_j p(t-jT_p).$$

Note that:

For each value of $j$, the pulse $a_j p(t-T_j)$ carries the information in its amplitude $a_j$. The amount of information depends on how many distinct values $a_j$ can take. Usually this number is represented by $M = 2^k$ for integer $k$, and we say that each pulse carries $k$ bits.
Since one pulse is transmitted every $T_p$ seconds, the pulse rate is $R_p = 1/T_p$ baud.

Now, (ignoring some more advanced considerations) the bandwidth of $s(t)$ is the same as the bandwidth of the pulse $p(t)$, denoted $B_p$. The pulse that has the minimum possible bandwidth for a given $R_p$ is a sinc pulse, and in that case the bandwidth of $s(t)$ is given by $$B_{min} = \frac{R_p}{2} = \frac{1}{2}\frac{1}{k}R_b.$$ Practical pulses achieve a fraction larger than $1/2$, or in other words, they require a larger bandwidth. If we define this fraction by $c$, we have $$B = c\frac{1}{k}R_b.$$

In conclusion, it seems to me that:

In your equations, $N$ is the bit rate, $r$ represents the number of bits transmitted per pulse, and $c$ represents the bandwidth efficiency of the chosen pulse. The variable name $k$ is more commonly used than $r$.
The first equation in your book is wrong; it should be $S = N \times \frac{1}{r}$.
The equation for the bandwidth is correct.


As stated in the paper citation, $(L-1)N$ zeros are padded in frequency domain to achieve an $L$ times oversampling, where $N$ is the number of subcarriers.

Example: Let $N=256$ and $L=4$. For the PAPR simulations use a 1024-IFFT. The 256 "center" inputs of the IFFT are the original OFDM frequency domain signal, the rest is zero. Thus $(L-1)N = 768\,$ zeros have been inserted.

"Nyquist-rate sampling" here means that all IDFT input bins are used and are in general unequal to zero, as for $L=1$ the number of zero padded subcarriers is 0. Since the IDFT input can be interpreted as a frequency domain signal, this means that the complete available bandwidth from $-f_\mathrm s/2$ up to $f_\mathrm s/2$ carries signal energy. In contrast, some of the available bandwidth is unused if zero padding in the frequency domain is applied. This is referred to as oversampling in analogy to sampling an analog signal at a higher rate than twice its maximum frequency. In the context of this question it is worth noting that the oversampling described in the article is only used for simulation purposes and that the signal will usually be analyzed at the IDFT output. A transmission of the oversampled signal is usually not done as the purpose of oversampling here is only PAPR analysis.