Group cascading IIR filters of different orders

You certainly can. In theory you can cascade any type of linear time invariant filter (FIR, IIR, overlap-add, warped allpass, etc.) of any filter order in any sequence and happily mix them up. The transfer function of the cascade is simply the product of all the individual transfer functions.

In practice, especially for fixed point systems, you need to watch gain, stability and noise, so depending on your application there may be some practical constraints.

The most common IIR implementation is cascaded biquad filters (also known as second order sections) since they are the easiest to manage in terms of noise, gain and stability.


What is the link between equation of a continuous signal versus equation of its sampled form?

The sequence $x[n]$ equals the continuous-time signal $x(t)$ sampled at $t=nT$, where $T=1/f_s$ is the sampling period:

$$x[n]=x(nT)=\cos(4000\pi nT)=\cos(4000\pi n/f_s)\tag{1}$$

So if $x[n]=\cos(n\pi/2)$ you just have to compare this argument to the argument in $(1)$ to figure out what $f_s$ is.

The normalized frequency of the discrete-time signal $x[n]$ is $\omega_0=\pi/2$ (in radians). It is related to the frequency $f$ in Hertz by

$$\omega_0=\frac{2\pi f}{f_s}$$


Discrete Cosine Transform (DCT) Coefficient Distribution

Probably the value of the DCT coefficient and its frequency (i.e., a histogram plot). However, it's impossible to say that with certainty given the information you have provided. If you found that plot in the linked paper, then the definitions of the x- and y-axes are in all likelihood given therein. The paper is behind a paywall, so it's not helpful to anyone without an IEEE account.
Yes. Check the scipy documentation for the dct() function.
See 1.
This is impossible to answer. What does appropriate mean in this context?
To classify a coefficient as either low, middle, or high frequency, compute its distance from the DC bin. For bin (m,n), its distance will be r = sqrt( (m-1)^2 + (n-1)^2 ). Judging by the picture you included, it looks like low frequencies are those with r < M/2, middle frequencies are those with M/2 <= r < M, and high frequencies are those with r >= M.

EDIT:

From the linked paper:

  Furthermore, it has been explained that noise addition to an image
  changes its DCT coefficient distribution [9], as we confirm in Fig. 2
  for a document image. Part (a) of the figure is an example of
  noiseless document image, part (b) is the noisy version of image (a)
  due to Gaussian blurring and Gaussian noise addition, and part (c)
  shows the higher frequency DCT coefficient distribution of of part(a)
  and (b) respectively. A lower-peaked histogram belongs to the
  ground-truth (GT) image, while a sharp-peaked histogram belongs to the
  noisy document image.

So part (c), which is the image you included, shows the histograms of the "higher frequency" DCT coefficients. Thus, the x-axis is the value of the DCT coefficient for the high frequency band, and the y-axis is its normalized histogram value. (The normalized histogram value is sometimes called "frequency", but that can get confusing in this case since DCT coefficients also represent frequencies. The normalized histogram value is just the number of occurrences of a certain value divided by the total number input samples to the histogram.)


Is sparse signal in time domain is sparse in Frequency domain

Your central question first:

  My question, does the results of FFT in frequency domain must be sparse since its equivalent time signal is sparse too ?

No. In fact, it does the opposite: signals sparse in one domain are typically dense in the other.

If you want so, that's happening for the same reason that we have Heisenberg's Uncertainty Principle: Operators like the DFT don't allow us to preserve correlation / concentration of energy when going from one domain to the other (or back).

Let me comment on your original statement:

  In OFDM system, if we have a multi-path environment, we supposed to have a sparse received signal in time domain, then FFT is performed to convert it into frequency domain.

That is false. The characterizing property of multipath environments is that they convolve the transmit signal with a nontrivial impulse response.

That converts even very sparse time-domain signals (e.g. a single dirac impulse) to dense ones. The only case where a channel would introduce sparsity would be if the the impulse response of the channel happened to be a matched filter to the transmit waveform.

But that can not be the case for an OFDM system, as that would directly contradict the principle of OFDM which is that you transport different bits of information on different frequency bins.

You might be considering a predistorting system, but then you'd mention that, and also, you'd know about the spectral properties of your signal (and hence could've negatively answered your original question yourself).

So, it feels like you've either taking something completely out of context, or you're misunderstanding something fundamentally.


What approximation techniques exist for the square super-root function?

Some numerical stabs in the dark yielded the following for an iterative approach:

We're looking for the solution y = f(x) where y^y = x.

$$y \ln y = \ln x$$

$$y = g(x,y) = e^{\frac{\ln x}{y}}$$

The value for $y$ is a fixed point of the above equation, and empirically it seems to converge for some values of $x$, but for larger values of $x$ it oscillates or diverges.

Then I tried an approach similar to Newton's iterative square-root:

$$y = \frac{y_{previous} + y_{*}}{2} = \frac{y + e^{\frac{\ln x}{y}}}{2}$$

where y* is supposed to represent a nonconverging but optimistic answer that maintains accuracy if you happen to guess an accurate initial value (in square root y2 = x, it's y* = x/y).

This appears to converge, but very slowly at the low end of $x$ (near $x_{min} = (\frac{1}{e})^{\frac{1}{e}}$ )

It also looks like a good initial guess is $y_0 = \ln(x)+1$.

So I figured maybe there's a better-converging solution:

$y = (1-a) \times y + a \times g(x,y)$ for some value of $a$ that is a function of $x$.

Then I found something interesting.

If I get a converging answer $y$ from the above approach for $y^y = x$, and then compute $y_2 = g(x,y+ \epsilon) = e^{\frac{\ln(x)}{y+\epsilon}}$, it appears as though $y_2-y$ = approximately $\epsilon \times (-\ln(y))$.... e.g. if we had a guess $y_1 = y + \epsilon$ for some unknown $\epsilon$, and computed $y_2 = g(x,y_1)$, then $(y_2 - y) \approx \epsilon \times (-\ln(y)) = (y_1 - y) \times (-\ln(y))$. (Just to clarify, I have no analysis to verify this, but the numbers just popped out of some numerical evaluation I performed.)

Solve for the linear terms in $y$, and you get $y = \dfrac{y_2 + \ln(y) \times y_1}{1 + \ln(y)}$... use $\ln(y_1)$ in place of $\ln(y)$ and you get this iterative approximation:

$$\begin{align*}
    y[n+1] &= \frac{g(x,y[n]) + \ln(y[n])\times y[n]}{1+\ln(y[n])}\\
           &= \frac{e^{\frac{\ln(x)}{y[n]}} + \ln(y[n])\times y[n]}{1+\ln(y[n])}
\end{align*}$$

This appears to work very well, with the initial guess $y = 1+\ln(x)$, and appears to converge within 4 or 5 iterations.

(Someone could probably show that this is equivalent to Newton-Raphson in some way, but I think it's beyond my capability.)


The Toeplitz matrix is used to compute correlation and convolution using matrix multiplication. Below is a graphic showing how to use a Toeplitz matrix specifically to perform convolution using matrix multiplication. A cross-correlation function can also be done following the same process by simply reversing the order of the coefficients in the multiplied vector (in this case we would multiply by [h3 h2 h1] to perform correlation).

This is clearer by observing the similarity and differences in the formulas for discrete time convolution and cross-correlation:

CONVOLUTION:

$$r[n]*h[h] = \sum_{m=-\infty}^\infty r[m]h[n-m]$$

CROSS CORRELATION:

$$\rho_{rh}[n] = \sum_{m=-\infty}^\infty r[m]h[m+n]$$

The code in the blue cell in the graphic is the implementation in either Matlab or Octave. This applies to auto-correlation or cross-correlation as an alternate to using the xcorr() function directly.

The utility of doing this is clearer in the following post, where I used such an operation to help demonstrate the underlying operation of the LMS equalizer:

Compensating Loudspeaker frequency response in an audio signal


What are the fundamental reasons that suggest us to minimise transmission power for energy-efficient radio link design?

Note that there are a host of practical secondary effects of increasing transmission power that do not scale with potentially less required transmission time, such as: battery impedance, voltage converter startup time, oscillator stabilization and carrier synchronization time, leakage or resistive losses due to larger required transistor sizes and/or higher voltage capacitors, more EMI shielding due to any higher digital edge rates required for processing, antenna shielding/spacing to limit human RF exposure, larger required heat spreaders/sinks, and etc.  The goal is often to minimize, not transmission power, but all these secondary power losses or inefficiencies, and/or their associated device costs, assuming some lower bandwidth or latency of delivering messages meets design goal targets.


Question on sketching the amplitude spectrum of signals

It's really very simple, but should have been clarified in your tutorial. It may turn out that a signal's Fourier Transform is real (I'll let you figure out when that happens). In that case, the spectrum $X(f)$ can be drawn in a single plot using positive and negative numbers, without having to split it up into magnitude and phase spectrums.

You can, of course, represent a real number in polar form. If the number is positive, its phase is zero, and if it is negative, its phase is $\pi$. What this means is that, whether you plot the real Fourier transform using negative and positive numbers, or you split it into magnitude and phase plots, you're conveying the exact same information.


How do you find the frequency and amplitude from a DFT?

For a lower quality frequency estimate, but still more precise than just using the (center) frequency of the peak magnitude DFT result bin, try using parabolic interpolation (using 3 neighboring points) of the spectral peak.

You can't get the phase from a DFT magnitude-only plot.  For phase you will need an estimator that uses the complex DFT results.

For a higher quality estimate, including a phase estimate, try using windowed-Sinc interpolation of the complex DFT result components, perhaps followed by successive approximation to find a better estimate of the magnitude peak near the peak DFT bin.  The quality of the estimate will vary with the width of the window.

Equivalent to Sinc interpolation would be to use a much longer zero-padded DFT or FFT to get more closely spaced frequency points.

There are a few added comments regarding the above frequency estimation methods on my blog:
http://www.nicholson.com/rhn/dsp.html#1

Note that near DC or Fs/2 (closer than the width of your interpolation kernel), you may need to use a mirrored interpolation kernel.  Note that hybrid methods are possible: e.g. using Sinc interpolation to up the resolution by 4X or 10X, then using parabolic or higher-degree polynomial interpolation to estimate the peak location, then using Sinc interpolation to estimate the magnitude and phase at that estimated peak frequency.  etc. Note that sinusoidal phase estimation in noise will often be more accurate when the phase is referenced to the center of the DFT window (N/2), not to sample point 0.  Referencing phase to the center of a data vector can be done using some form of fftShift.

The question does not ask about frequency estimators that do not involve a DFT, but note that in zero noise, only 3 or 4 unaliased points are required to solve for 3 unknowns (frequency, magnitude and phase of a pure unmodulated sinusoid).  See: http://www.claysturner.com/dsp/3pointfrequency.pdf and http://www.claysturner.com/dsp/4pointfrequency.pdf


Will the capacity of a channel becomes unbounded if i increase its signal-to-noise ratio $S/N$ without limit?

By the sampling theorem, uniform samples taken at more than twice the bandwidth of a band-limited signal can be used to reconstruct the signal perfectly by sinc interpolation. For each different set of sample values there is a different band-limited function. In absence of noise, by increasing the bit depth of the samples one can store an unlimited amount of information per sample into the band-limited signal, and get it back simply by sampling and digitizing it again.


OFDM time vs. frequency domain channel estimation/equalization

Traditionally, OFDM became popular in WiFi and LTE because the channel model consisted of multi-path. That is, the radio signal transmitted in 1-6GHz frequencies bounced from various obstacles (walls, trees, cars, humans) at the receiver. Of course this is time varying because obstacle position or transmitter/receiver position also changes. But to simplify the computation, if we consider a multi-path channel at time $t$, its baseband model for $L$ tap multipath is

$$
h[n] = \sum_{k=0}^{L-1}a_k\delta[n-k]
$$
where $a_k = r_ke^{j\phi_k}$ is the complex tap for path $k$. If you plot the FFT of this channel, you will find it is not a flat but rather having different gains at different frequencies. To equalize this in time domain, you would use adaptive filtering techniques. The problem is complexity and time taken for adaptive filtering to estimate the channel. By the time adaptive filter error tapers down, channel would have changed (due to movement of transmitter/receiver/obstacle). So for each OFDM burst, there is a need to immediately equalize to demodulate the packet and ultimately show the data on the device! This is where advantage of OFDM comes into picture. Even though the complexity of OFDM is higher (overhead of preamble, pilot symbols, need to maintain orthogonality), the equalization simplifies to a single tap equalizer in frequency domain as the whole frequency selective channel is split into small flat-fading channel due to IFFT -> Cyclic Prefix -> FFT method. You make sure the circular convolution happens by inserting cyclic prefix so in frequency domain it turns into a point-by-point equalization $Y[k]/H_{eq}[k]$ where $Y[k]$ is the received symbol at subcarrier $k$, and $H_{eq}[k]$ is the equalizer tap you computed using pilots.

EDIT: I understand the channel model given above is very simplified view of multi-path model. The simplification is explained in the famous book by Prof. David Tse and Prof. Pramod V (https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_chapter2.pdf)


Relationship Between Sampled Continuous and Discrete Time Signals

You've shown that the (continuous-time) Fourier transform (CTFT) of a sampled continuous-time signal equals the discrete-time Fourier transform (DTFT) of the corresponding discrete-time signal. In both cases, the spectrum is periodic. The independent variables are related by

$$\omega=\Omega T\tag{1}$$

where $\Omega$ is the angular frequency ($\Omega=2\pi f$), $T=1/f_s$ is the sampling period ($f_s$ is the sampling frequency), and $\omega$ is the normalized angular frequency, which is often used as the independent frequency variable when describing spectra of discrete-time signals. Note that you could also use $\Omega$ as the independent variable for discrete-time signals, but it always appears together with the sampling interval $T$. This is an important property of sampled or discrete-time signals, and that's also why the two variables $\Omega$ and $T$ are commonly combined into a single variable $\omega=\Omega T=2\pi f/f_s$.

Note that sampling results in a periodization of the spectrum, and this is why information can be lost. This is true for both spectra $X_s(j\Omega)$ and $X(e^{j\omega})$. The periodicity of the spectra implies that the fundamental frequency interval $f\in [-f_s/2,f_s/2]$, or, equivalently, $\Omega\in [-\pi/T,\pi/T]$ or $\omega\in [-\pi,\pi]$ contains all information of the sampled signal, and all frequencies outside that interval are redundant. This is an inherent property of sampled (or  discrete-time) signals, and it is the cause of aliasing, i.e. the mapping of different continuous-time signals to the same sampled signal.


improving pre-processing (noise removal) orange seven segment display for OCR reading

The shapes you are interested in are very regular while the noise you are trying to deal with will be less regular. Ideally, you will only have two shapes at the end of your processing: wide rectangles for the horizontal parts of the display and tall rectangles for the vertical parts. You can leverage this fact in how you shape your structuring element when you use opening to remove noise.

This code is for Matlab (or Octave, its free cousin), but should easy to implement in Python. Note that I screen grabbed the image from your GitLab file and resized it - so these dimensions likely won't match yours exactly. I also inverted the image so the foreground is white. Also, the images below were resized to 400% before posting.

pkg load image
img = imread("download.png");
img = img * 255;
imshow(img)

First, use a wide structuring element for opening so only wide rectangles are kept.

SE_h = strel('rectangle',[1,8]);
h_lines = imerode(img, SE_h);
h_lines_reconstructed = imdilate(h_lines, SE_h);
imshow(h_lines_reconstructed)

The repeat for tall rectangles:

SE_v = strel('rectangle',[8,1]);
v_lines = imerode(img, SE_v);
v_lines_reconstructed = imdilate(v_lines, SE_v);
imshow(v_lines_reconstructed)

Recombine and display:

combined = or(v_lines_reconstructed, h_lines_reconstructed);
imshow(combined);


Read the audio data

[y,Fs,bits] = wavread('file1.wav');

Apply window function on that data

ham=hamming(Nsamps);
windowed = y .* ham;

Apply fft on the result of window function

ham_fft = fft(windowed);
ham_fft = ham_fft(1:Nsamps/2);      %Discard Half of Points

Make power spectrum of fft

PowSpec = abs(fft(ham_fft)).^2;

divide power spectrum in Barks bands

So here maybe start interesting things, based in wikipedia the 24 Critical Barks Bands in Hertz is:

20,100,200,300,400,510,630,770,920,1080,1270,1480,1720,2000,2320,2700,3150,3700,4400,5300 6400,7700,9500,12000,15500

Now to split the power spectrum in Barks bands you need to do something like:

BandBarks = [20 100 200 300 400 510 630 770 920 1080 1270 1480 1720 2000 2320 2700 3150 3700 4400 5300 6400 7700 9500 12000 15500];

for l = 1:length(BandBarks)-1
    CurrentBand = [BandBarks(l):BandBarks(l+1)];
    PowSpec(CurrentBand,:) %% here is the power spectrum for each Barks bands
end

OK by your commented error seems that your audios are in stereo, convert to mono, or get just one channel!

my example above assumes that each BIN number is the equivalent frequency (course this is nothing common was just an example), to find what bin is your frequency you need calculate:

$$ BandBarks * (\frac{Fs}{length(FFT)})$$

So based in this equation:

Bins = zeros(length(BandBarks),1);
for j = 1:length(BandBarks)

    Bins(j)=round(BandBarks(j)/(Fs/length(FFT)))

end

Now here how to plot individually each Bark band:

for l = 1:length(Bins)-1
    CurrentBand = [Bins(l):Bins(l+1)];
    figure(l);
    plot(PowSpec(CurrentBand,:));
end

I think I can not be clearer than that


Is the inverse of a causal system also causal?

It's not sufficient to only consider causality, you also need to check whether the inverse system is stable, otherwise it can't be implemented. If $G(z)$ has zeros on the unit circle, it cannot be inverted. If $G(z)$ has no zeros on the unit circle, but if there are zeros outside the unit circle, then there is no causal and stable inverse, because the zeros of $G(z)$ become the poles of the inverse system. There is, however, a stable inverse but it is not causal. If $G(z)$ has all its zeros (and poles) inside the unit circle, then there exists a causal and stable inverse (with all its poles and zeros inside the unit circle). Such a system is called strictly minimum-phase, and its inverse is also a strictly minimum-phase system.


no.  but, for a causal LTI digital filter with a rational transfer function, it is necessary that the number of zeros, $M$, not exceed the number of poles, $N$.

$$\begin{align}
\\
 H(z) \triangleq \frac{Y(z)}{X(z)} & = H_0 \frac{(z-q_1)(z-q_2)(z-q_3)...(z-q_M)}{(z-p_1)(z-p_2)(z-p_3)...(z-p_N)} \\
\\
 & = H_0 \frac{z^M (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})}{z^N (1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = z^{M-N} \ \frac{H_0 \ (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})}{(1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = z^{M-N} \ \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_M z^{-M} }{ 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... +  + a_N z^{-N} }
\end{align}$$

this gets manipulated:

$$\begin{align} Y(z) & \left(  1 + a_1 z^{-1} + a_2 z^{-2} + ... + a_N z^{-N}  \right) \\
& \quad = X(z) z^{M-N} \left( b_0 + b_1 z^{-1} + b_2 z^{-2} + ... + b_M z^{-M} \right) \end{align}$$

or

$$\begin{align} Y(z) = \ & b_0 z^{M-N} X(z) + b_1 z^{M-N-1} X(z) + b_2 z^{M-N-2} X(z) + ...  \\
 & \quad ... + b_M z^{-N} X(z) - a_1 z^{-1}Y(z) - a_2 z^{-2}Y(z) - ... - a_N z^{-N}Y(z) \end{align}$$

inverting the Z Transform, this becomes

$$\begin{align} y[n] = \ & b_0 \ x[n+M-N] + b_1 \ x[n+M-N-1] + b_2 \ x[n+M-N-2] + ...  \\
 & \quad ... + b_M \ x[n-N] - a_1 \ y[n-1] - a_2 \ y[n-2] - ... - a_N \ y[n-N] \end{align}$$

for a causal system or filter, we cannot look into the future for any input sample $x[n]$.  that means $M \le N$ so that there is no negative delay on the $x[n+M-N]$ term.  a negative delay means the LTI filter is peeking into the future, which is not causal.

but if you remove that leading delay factor and have this:

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_M z^{-M} }{ 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... +  + a_N z^{-N} } $$

there is no reason that $M$ cannot be larger than $N$. but the number of poles will still be as large as $M$. if $M>N$, then some of the poles will have to be at the origin $z=0$. in fact, if $N=0$, then you have an FIR filter and all of the poles are at the origin.


Blind source separation on images - with a known source

You already know one of the sources here, so you don't need ICA.

There are infinitely many couples $k$ and $B$ that verify:
$$
A = B + kC,
$$
 but in your case if you assume that $B$ and $C$ are uncorrelated, we can choose the solution that minimizes the correlation between $B$ and $C$.

Using the previous equation and the linearity of the correlation operator, we have:
$$
\text{corr}(A,C) = \text{corr}(B,C) + k \text{corr}(C,C)
$$

So you need to find $k$ that minimizes $(\text{corr}(A,C) - k \text{corr}(C,C))^2$, which is simply:

$$
k = \frac{\text{corr}(A,C)}{\text{corr}(C,C)}.
$$
And finally:
$$
B = A- \frac{\text{corr}(A,C)}{\text{corr}(C,C)}C
$$


The "poles-inside-unit-circle" stability criterion only applies to causal systems. Your system is not causal because it uses one sample from the future owing to the $z$ term.

The general technique to check for stability involves looking at the regions of convergence (ROC) of $H(z)$. If the ROC includes the unit circle, then the system is stable.

See also this answer to the question: What is the relationship between poles and system stability?


How does one compute the band width of a signal in the time domain $x(t)$?

As far as I know there is no way to get $f_{max}$ and $f_{min}$ directly from time domain (zero crossing rate will not help if the signal has a non-zero DC offset for example) you need to look at the frequency domain.

If there is just one term, the smaller frequency is just $0$

Given a time signal and its sampling rate ($Fs$), you can take the discrete Fourier Transform to study its frequency spectrum. Unfortunately you cannot use this to decide your sampling rate if your initial sampling rate ($Fs$) is lesser than the 2$*$bandwidth (2$*$[highest frequency] if you are looking for perfect reconstruction) as you have already lost the data during initial sampling(look up aliasing for further reading).

Bottom Line: You need to have an estimate in what frequency range your signal will lie so that you could sample it accordingly.


Difference between Sub Sampling and Down Scaling of Images

Have you looked to the pyUp and pyDown functions in opencv? They upscale and downscale images using a gaussian pyramids.

opencv image filtering page


The filter is a linear time-invariant (LTI) system that adds memory to the model. Combined with a memoryless nonlinearity, the filter allows to model non-linear systems with memory, unlike a pure AM-AM and AM-PM model, which models a memoryless nonlinearity. An ARMA-filter is the same as an infinite-impulse response (IIR) filter, which computes its output as a linear combination of weighted and delayed versions of the input and output signal. IIR filters need to be implemented by recursive structures.

An ARMA (or IIR) filter is a discrete-time (linear time-invariant) filter whose response $y[n]$ to an input sequence $x[n]$ is computed as follows:

$$y[n]=\sum_{m=0}^Mb_mx[n-m]+\sum_{l=1}^Na_ly[n-l]\tag{1}$$

where $b_m$ and $a_l$ are the filter coefficients, which are chosen such that a certain desired filter characteristic (e.g., low pass, high pass, etc.) is achieved.


Consider you have a bandwidth of B. The sampling rate $F_s$ which is related to B using the relation $F_s > 2B$. Now the sample time is $\frac{1}{F_s}$. So if you have N samples  then total time is $\frac{N}{F_s}$

Is that what you are looking for?


How to turn Morlet wavelet into second generation wavelet?

The Lifting scheme is special algorithm that allows a very efficient computation of the discrete wavelet transform (DWT), via special filters (QMF) and down/up-samplers. In order to compute a DWT the wavelet has to fulfill orthogonality constraints under translation and dilation.

The Morlet wavelet does not fulfill these contraints. Thus, you cannot compute a DWT with it und consequently you cannot use the lifting scheme. The Morlet wavelet is only applicable for computing a continuous wavelet transform. If you need a fast approximation you can have a look here: Fast Gabor-Like Windowed Fourier and Continuous Wavelet Transforms


Confused about the Fourier Transform of $e^{at}u(t)$

$\omega$'s value does nothing for the convergence of the integral. $a$'s does.

$a$ is not complex, it's real only.

If $a$ is 0, the integral doesn't converge.


Neither. You calculate the power spectrum from the signal directly. Both RMS value and peak are a single value that quantifies a long period. Hence, for all frequencies higher than the reciproce of that period will be lost.


correlation matrix for complex white noise does not equal to identity matrix?

Why don't you simply compute $X^{H}  X / N$ ?
clear all;
close all;

fs = 4000;
N = 100*fs;
R = 0.5;% overlap percent
Nseg = fs;

M = 8;

I noticed your data was not complex, since you metioned complex
in the question here I made it complex.
X = 0.5*randn(N,M) + 0.5j *randn(N,M);

This will be close to the identity matrix there will be small complex perturbations
the expected perturbations are inversely proportional to $\sqrt{N}$. The $X^H$ corresponds to X', it is the conjugate transpose of $X$.
R_hat = X' * X / N

Edit
Here I compute RR such that RR(:,:,k) is the accumulated covariance matrix of the frequency bin $k$ for the different inputs.
clear all;
close all;

fs = 4000;
N = 100*fs;
R = 0.5;% overlap percent
Nseg = 256;

M = 8;

for m = 1:M
   x(:,m) = randn(N,1) ;
end

segNum = floor((N-R*Nseg)/((1-R)*Nseg));

RR = zeros(M, M, Nseg);

for n = 1:segNum
   y = x([1:Nseg]+(n-1)*Nseg*(1-R),:) ;
   Y = fft(y);
   for k = 1:Nseg
     %% the correlation for for one bin one slice.
     RR(:,:,k) = RR(:,:,k) + Y(k, :)' * Y(k, :);
   end
end
%% Normalize
RR = RR * (1/segNum/Nseg);

%% Acculate over all bins
R_hat = 0;
for n = 1:Nseg
   R_hat = R_hat + RR(:,:,k);
end

R_hat = R_hat/Nseg


How to produce uniformly distributed low-frequency noise?

Well, this is one way you could do it: with LTspice, use the white() function. There are also rand() ([0,1] V, pulses) and random() ([0,1] V, smooth pulses), but white(2*fs*time) will give you what you want between [-0.5, 0.5] V:

B1 is a behavioural current source (bi, or bi2, as it appears in the component selection dialog, F2), and the .wavecommand exports it as a .wav file named output.wav, with 16 bits resolution (can be anything from 1 to 32) and 44.1 kHz sampling frequency (can be any number from 1 to 4294967295). For this case, the file is 8,820,046 bytes, which can be imported. This is a 2 Hz white noise, band limited, properly sampled signal, 100 sec long. Optional settings can include a tighter timestep, though opt plotwinsize=0, here, disables waveform compression and for a 2 Hz signal, it may not make much sense.


Subtracting signals to re-construct the middle

No, you can't restore a low - frequency signal by DFT of its parts due to the periodic DFT nature and because
$$
IDFT(DFT( \frac{U(t) + V(t)}{2})) == \frac{U(t) + V(t)}{2}
$$
by definition.
You can apply filtering or smoothing to remove the high - frequency noise.


Perfect reconstruction: what properties must the analysis LPF satisfy

So the central question here is what the restrictions for the filters of this class of filterbanks (i.e. $N=2$-filter bank with the HPF generated from a lowpass/highpass transform of the LPF) is.

Let us for that first introduce a couple of helping definitions:

$$\begin{align}
\mathbf{h}(z) &= \begin{bmatrix}h_0(z)\\h_1(z)\end{bmatrix}&\text{analysis filterbank matrix}\\
\mathbf{r}(z) &= \begin{bmatrix}r_0(z)\\r_1(z)\end{bmatrix}&\text{reconstruction filterbank matrix}\\
\end{align}$$

With these in place: Let's consider the polyphase TypeI representation $A(z)$ (analysis) and TypeII representation $R(z)$ (reconstruction) of these filter banks (by the way, this is why we have these different types of polyphase deconstructions); I hope my notation is clear, $h_{1,0}$ is the zeroth polyphase component of the high pass filter.

$$\begin{align}
\mathbf{h}(z) &=
   \begin{bmatrix}
     h_{0,0}^{\mathrm I}(z^2) & h_{0,1}^{\mathrm I}(z^2)\\
     h_{1,0}^{\mathrm I}(z^2) & h_{1,1}^{\mathrm I}(z^2)
   \end{bmatrix}
   \begin{bmatrix}
     z^{-0}\\
     z^{-1}
   \end{bmatrix}\\
&= \mathbf A(z^2) \mathbf e_z(z)\\
\mathbf{r}^T(z) &=
   \begin{bmatrix}
     z^{-1} &
     z^{-0}
   \end{bmatrix}
   \begin{bmatrix}
     r_{0,0}^{\mathrm {II}}(z^2) & r_{0,1}^{\mathrm {II}}(z^2)\\
     r_{1,0}^{\mathrm {II}}(z^2) & r_{1,1}^{\mathrm {II}}(z^2)
   \end{bmatrix}\\
&= \mathbf{ \tilde e}_z(z) \mathbf R(z^2) \\
\end{align}$$

Now, looking at this, we simply define perfect reconstruction "intuitively" as the case where applying $\mathbf R$ to the output of $\mathbf A$ leads to the input of $\mathbf A$, of course delayed, and maybe scaled with a constant $c$:

$$\begin{align}
  h,\,r &\text{ form a perfect reconstruction fb}\\
  &\iff\\
  \mathbf R(z)\mathbf A(z) &= c z^{-k} \mathbf I_N & c \text{ const.},\, k\in\mathbb N
\end{align}$$

That gives us a clear method to move forward:

$$\begin{align}
  \mathbf R(z)\mathbf A(z) &=    \begin{bmatrix}
     h_{0,0}^{\mathrm {II}}(z) & h_{0,1}^{\mathrm {II}}(z)\\
     h_{1,0}^{\mathrm {II}}(z) & h_{1,1}^{\mathrm {II}}(z)
   \end{bmatrix}
   \begin{bmatrix}
     r_{0,0}^{\mathrm I}(z) & r_{0,1}^{\mathrm I}(z)\\
     r_{1,0}^{\mathrm I}(z) & r_{1,1}^{\mathrm I}(z)
   \end{bmatrix}\\
&= \begin{bmatrix}
h_{0,0}^{\mathrm {II}}(z)r_{0,0}^{\mathrm I}(z) +  h_{0,1}^{\mathrm {II}}(z)r_{1,0}^{\mathrm I}(z) &
h_{0,0}^{\mathrm {II}}(z)r_{0,1}^{\mathrm I}(z) +
h_{0,1}^{\mathrm {II}}(z)r_{1,1}^{\mathrm I}(z)\\
h_{1,0}^{\mathrm {II}}(z)r_{0,0}^{\mathrm I}(z) +
h_{1,1}^{\mathrm {II}}(z)r_{1,0}^{\mathrm I}(z) &
h_{1,0}^{\mathrm {II}}(z)r_{0,1}^{\mathrm I}(z) +
h_{1,1}^{\mathrm {II}}(z)r_{1,1}^{\mathrm I}(z)
\end{bmatrix}\\
&\overset!= \begin{bmatrix}c z^{-k}&0\\0&c z^{-k}\end{bmatrix}
\end{align}$$

EDIT now, let's do our jobs of inserting our conditions:

$$r_0  = h_0 (L-1-n)$$

notice that the substitution of $n$ by $L-1-n$ implies that the typeI polyphase decomposition of $r_0$ is the same as the typeII decomposition of $h_0$!

Hence, the diagonal entries of above matrix get a lot easier to calculate:

$$\begin{align}
  \mathbf R(z)\mathbf A(z) &=\begin{bmatrix}
h_{0,0}^{\mathrm {II}}(z)r_{0,0}^{\mathrm I}(z) +  h_{0,1}^{\mathrm {II}}(z)r_{1,0}^{\mathrm I}(z) &
h_{0,0}^{\mathrm {II}}(z)r_{0,1}^{\mathrm I}(z) +
h_{0,1}^{\mathrm {II}}(z)r_{1,1}^{\mathrm I}(z)\\
h_{1,0}^{\mathrm {II}}(z)r_{0,0}^{\mathrm I}(z) +
h_{1,1}^{\mathrm {II}}(z)r_{1,0}^{\mathrm I}(z) &
h_{1,0}^{\mathrm {II}}(z)r_{0,1}^{\mathrm I}(z) +
h_{1,1}^{\mathrm {II}}(z)r_{1,1}^{\mathrm I}(z)
\end{bmatrix}\\
&=\begin{bmatrix}
h_{0,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm {II}}(z) +
h_{0,1}^{\mathrm {II}}(z)-h_{1,0}^{\mathrm{II}}(z) &
h_{0,0}^{\mathrm {II}}(z)h_{0,1}^{\mathrm {II}} +
h_{0,1}^{\mathrm {II}}(z)-h_{1,1}^{\mathrm{II}}(z)\\
h_{1,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm{II}}(z) +
h_{1,1}^{\mathrm {II}}(z)-h_{1,0}^{\mathrm{II}}(z) &
h_{1,0}^{\mathrm {II}}(z)h_{0,1}^{\mathrm{II}}(z) +
h_{1,1}^{\mathrm {II}}(z)-h_{1,1}^{\mathrm{II}}(z)
\end{bmatrix}\\
&=\begin{bmatrix}
h_{0,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm {II}}(z) -
h_{0,1}^{\mathrm {II}}(z)h_{1,0}^{\mathrm{II}}(z) &
h_{0,0}^{\mathrm {II}}(z)h_{0,1}^{\mathrm {II}} -
h_{0,1}^{\mathrm {II}}(z)h_{1,1}^{\mathrm{II}}(z)\\
h_{1,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm{II}}(z) -
h_{1,1}^{\mathrm {II}}(z)h_{1,0}^{\mathrm{II}}(z) &
h_{1,0}^{\mathrm {II}}(z)h_{0,1}^{\mathrm{II}}(z)
h_{1,1}^{\mathrm {II}}(z)h_{1,1}^{\mathrm{II}}(z)
\end{bmatrix}\\
\end{align}$$

To achieve the zero entries in the identity matrix, it is necessary that

$$\begin{align}
h_{0,1}^{\mathrm {II}} h_{0,0}^{\mathrm {II}}(z)&= h_{0,1}^{\mathrm {II}}(z)h_{1,1}^{\mathrm{II}}(z)\\
h_{1,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm{II}}(z) &= h_{1,0}^{\mathrm{II}}(z)h_{1,1}^{\mathrm {II}}(z)
\end{align}$$

This breaks down to filters that need to have zeros at every $N$th position, since the 2-polyphased decomposition of a HPF/LPF transformed HPF is either the same as the original's decomposition, or the original's decomposition with a factor of $-1$.


How to sketch this multirate system spectrum?

You are on the right track.

  After upsampling the spectrum of X would look the same except the x axis would be divided by 2. So the limit of the spectrum would be at pi/4.

Yes, that's correct.

  Since the filter has amplitude 2 would the amplitude of X be increased from 1 to 2 or will the amplitude of X remain unchanged?

Recall that in the frequency domain, the application of a linear system transfer function $H_B(\Omega)$ is done via straight multiplication of the input signal's spectrum. So yes, since the filter has a magnitude of 2, the amplitude of its output is doubled.

  Finally, for the multiplication of the cosine it would be equivalent to multiply with 2 deltas located at -pi/2 and pi/2. Would the multiplication with this deltas cause a kind of mirroring effect on the spectrum?

Here's where I will make a small correction: multiplication by a cosine in the time domain yields convolution with a pair of Dirac impulse functionss in the frequency domain. Convolution with a Dirac impulse in the frequency domain is equivalent to a shift along the frequency axis. So, you end up with 2 copies of the filter's output spectrum, centered at $\omega = \pm \frac{\pi}{2}$.


Does sampling ~> Nyquist frequency capture harmonics?

This is called "aliasing" and it is simply a result of periodicity of sine.

$\sin ak$ = $\sin (a+2b\pi)k$ (for integers $k$,$b$ - $k$ being the sample index)

So, if you maybe have $a=1.3\pi$, you could also see it as a signal at $-0.7\pi$ (a signal of negative phase at $0.7\pi$), or rather: you will not see any difference in the samples.


Can one generate a signal with an asymmetric power spectral density?

@user129412, there are many ways to generate both continuous signals and discrete sequence having asymmetric PSDs. One way, for discrete sequences, is shown in the following diagram.

$x_r(n)$ is a narrow bandpass signal having spectral energy in the vicinity of $±f_c$ Hz. The complex-valued $x_c(n)$ output sequence is your desired asymmetrical sequence having spectral energy only in the vicinity of $+f_c$ Hz. (If your Hilbert transformer has M taps, where M is odd, then the upper path delay must be (M-1)/2 samples to achieve time synchronization between the two signal paths.)  If you desire sequences whose asymmetrical spectra are centered at zero Hz, many ways to obtain such sequences can be found at:

https://www.dsprelated.com/showarticle/153.php


Does Wav format supports writing complex values

You could create a standard 2 channel wav file. 1 channel for the real part and 1 for the imag. Wav files supports both 16bit integer, 32bit integer and 32bit floats etc. The channel data will be interleaved (I believed it's called). Meaning that your data will be written: [REAL][IMAG][REAL][IMAGE]...
Prolly just like your own bin file.

Standard programs like Audacity will be able to read the file and make calculations on each (or both) channels.


Help in designing filter for DUC (Digital Up-Converter) and DDC (Digital Down-Converter) using the built-in DUC and DDC blocks in MATLAB Simulink

Help design my system would be far more detail than what we typically provide here, but I can give you some suggestions that may help or at least lead to a more concise question.

  I do not get correct BER, which means my signal is getting aliased as
  I receive through the DDC.

Not getting the correct BER does not necessarily conclude this.

In your case the error rate appears to be 50% which suggests your receive data is completely independent of your transmitted data - one possible reason this is if the delay for all the filtering hasn't been properly accounted for when comparing the transmit to the receive signals. Instead of comparing the constellations, consider comparing the baseband modulated signals in the time domain as waveform amplitude versus time to help determine what the actual failure mechanism is. That will show you if everything is fine and you just have the expected time offset, or otherwise give you further visibility into where exactly the problem is occurring (look at the in/out of each block in the time domain).

Once you isolate the block that is not performing as you would expect it to, it will either be quite clear to you what the problem is or you will have a more concise question to ask here.

Good luck!


What would be the variance for complex number?

I will focus on the reason of the factor $1/2$ and leave aside the estimation things.

The exact understanding should be : if a scalar Gaussian random variable (rv) is circular symmetric, its real and imaginary parts must be uncorrelated (this is equivalent to independence if they are assumed jointly Gaussian) and identically distributed with zero mean. Thus, your Matlab code is correct for a rv $\sim \mathcal{CN}(0,1)$.

The story behind is that a complex random variable (rv) is simply a vector of two real random variables. A vector of $n$ complex rv is indeed a vector of $2n$ real rv.

You are talking about the case $n=1$, 1 complex Gaussian rv $Z = Z_r + jZ_i$ or a vector of 2 real Gaussian rv $[Z_r,Z_i]^T$. As with real Gaussian rv which is described by its variance, a real vector $[Z_r,Z_i]^T$ must be described by its covariance matrix
$$\mathbb{E}\left\lbrace [Z_r,Z_i]^T \times [Z_r,Z_i] \right\rbrace = \mathbb{E}\left\lbrace\begin{bmatrix}
        Z_r^2 & Z_rZ_i \\
        Z_iZ_r & Z_i^2 \\
        \end{bmatrix}\right\rbrace$$

Take a look at the variance $\mathbb{E}\left\lbrace ZZ^H \right\rbrace = \mathbb{E}\left\lbrace Z_r^2+Z_i^2 \right\rbrace$ and pseudo-variance $\mathbb{E}\left\lbrace ZZ^T \right\rbrace = \mathbb{E}\left\lbrace Z_r^2-Z_i^2+j2Z_rZ_i \right\rbrace$, the covariance matrix of real vector (or complex scalar) is fully described by the variance and pseudo-variance. Threfore, you need both variance and pseudo-variance to characterize a complex Gaussian rv (with the prior condition that its real and imaginary parts are jointly Gaussian).

Now we use the circular symmetry property : $e^{j\phi}Z$ has the same probability distribution as $Z$ for all real $\phi$. This leads to $\mathbb{E}\left\lbrace e^{j\phi}Z(e^{j\phi}Z)^T \right\rbrace = \mathbb{E}\left\lbrace e^{j2\phi}ZZ^T \right\rbrace = \mathbb{E}\left\lbrace ZZ^T \right\rbrace$ for all $\phi$ thus $E[ZZ^T] = 0$ and variance $\mathbb{E}\left\lbrace ZZ^H \right\rbrace$ is sufficient statistic for $Z$. Note that $\mathbb{E}\left\lbrace ZZ^T \right\rbrace = \mathbb{E}\left\lbrace Z_r^2-Z_i^2+j2Z_rZ_i \right\rbrace \implies \mathbb{E}\left\lbrace Z_r^2 \right\rbrace = \mathbb{E}\left\lbrace Z_i^2 \right\rbrace$ and $\mathbb{E}\left\lbrace Z_rZ_i \right\rbrace=0$, the real and imaginary parts are uncorrelated then independent (because they are Gaussian), with the same variance, this is the reason of the factor $1/2$.

To sum up, your code is correct because you are estimating circular symmetric complex Gaussian rv. The jointly Gaussian assumption between real and imaginary parts must be used. If this is not about circular symmetric rv (or real random vector with two elements), you must calculate the pseudo-variance.

For more details and to understand the formula of the wikipedia article, you can read this article Circular Symmetric Gaussian R.Gallager.


Averaging is filtering with a rectangular filter kernel, which has a poor frequency response, both in pass band ripple and stop band attenuation, Because of the poor stop band response it won't anti-alias very well.  The non-flat passband will distort the final result even if there is no high-frequency spectral content to alias.


What are Chunks when recording a voice signal

So essentially, @porten pointed out what chunks of sound are.

Let's look at your code:

CHUNK = 1024
# What is CHUNKS here ?

The chunk is like a buffer, so therefore each buffer will contain 1024 samples, which you can then either keep or throw away. We use CHUNKS of data, instead of a continuous amount of audio because of processing power. (Let's assume this was a continuous flow of data being received from a Microphone and is being recorded and saved) then it would just eat up the processor, thus causing potential crashes. In terms of Raspberry Pi / Arduino development (Where the RAM is very small) CHUNKING the data like this makes the stream flow easier and thus prevents memory leaks.

Another reason:

Let's assume that you wanted to implement an algorithm for determining whether something is speech, or, is just noise. How would/could you do this using a constant stream flow of sound data? It would be very difficult. Therefore, by storing this into an array (or list in your case) you can perform analysis on this data, for example RMS. Then you could have some threshold to determine if you want to keep the data, or the data is no good. E.g.

CHUNK = 1024

var input[CHUNK] = readAudio(CHUNK); // reads in 1024 of data

// compute RMS
a = RMS(input);
if(a >= THRES_VALUE)
{
    actual += input;
}

So in this example above, we record 1024 bytes of data, we then compute RMS and check to see whether or not the sound is speech or just white-noise.

To confirm as well, the chunks contain the data.


Why is the ideal low pass filter not achievable?

General Background: Representing a signal numerically in a computer requires the signal to be discrete. That means it's periodic in the other domain. When you use something like a DFT you need the signal to be discrete in both domains: That means the signals are ALSO periodic in both domains.

If I apply an ideal low pass filter, I should get an infinite sequence of data points because the sinc function ...

If you want to do this numerically in a computer, it's NOT a sinc function. The lowpass filter MUST be periodic in both domains and so what you get in the time domain is a periodic repetition of a sinc function with a time domain aliasing.
It's still infinite it time. In fact it MUST be infinite in time since it's periodic. As long as all signals involved are actually periodic with a shared Period, than things work "as expected" .
At this point the term "low pass" filter becomes complicated. Since the spectrum is periodic, it will pass arbitrarily high frequencies. So it needs to be interpreted: it removes higher frequencies in the base band or something like this.

$y\left(m\right)\ =\ \displaystyle\sum_{i=1}^{N}\cos\left(m\alpha_{i}\right)$

Your missing phase or sine term here but you are re-inventing the DFT or cosine transform. You run into exactly the same problem. How to choose N ?

Also if I knew that the signal had exactly N points

That's a contradiction. The signal that you described in your original equation is infinitely long. As soon as you make it finite in time, you are windowing it and that changes the spectrum and the math behind considerably.

then the signal could be represented at most by N real frequencies

Correct. That's for example the Discrete Cosine Transform, which is a close cousin of the DFT. But it makes no material difference here.

Then why is the ideal low pass filter not realizable?

It is realizable as long as everything periodic "the right way". The periodicity helps dealing with the infinite length in time: you need only to compute what's happening in one period and you have all information that can be had since all other periods are the same.
However most real world signals are NOT periodic "the right way", so ideal low pass filtering is rarely feasible.


Your diagram is incorrect. Computing the 2D-FFT of a 8x8 matrix is not the same as computing the 1D-FFT of a 1x64 vector. A 2D-FFT is "probing" your signal with templates which have both an horizontal and vertical frequency. A signal with a low horizontal frequency and a null vertical frequency is a progressive horizontal gradient. A signal whose vertical and horizontal frequency are both the nyquist frequency is a checkerboard pattern...

If you take an image and perform a 2D-FFT on it, you will find in the resulting matrix the following elements:

Cell (1, 1): DC component (average of all values in the matrix).
Cell (2, 1): first horizontal harmonic; frequency 1/8.
Cell (3, 1): second horizontal harmonic; frequency 2/8.
Cell (4, 1): third horizontal harmonic; frequency 3/8.
Cell (5, 1): fourth horizontal harmonic; frequency 4/8 (Nyquist).
Cell (6, 1): third horizontal negative freq. harmonic; frequency -3/8.
Cell (7, 1): second horizontal negative freq. harmonic; frequency -2/8.
Cell (8, 1): first horizontal negative freq. harmonic; frequency -1/8.
Cell (1, 2): first vertical harmonic; frequency 1/8.
Cell (1, 3): second vertical harmonic; frequency 2/8.
...
Cell (5, 5): fourth horizontal/vertical harmonic (Checkerboard pattern).
...

A good way to visualize this is to compute the 2D IFFT of a single pixel at each of the 64 positions of the matrix, and look at the real part.


Help with algorithm for modulating oscillator pitch using LFO

As said by Jason, this might simply be that you are not implementing your oscillator correctly - for example by multiplying the frequency by time instead of integrating it.

Note also - and this is unrelated to your topic but really worth observing - that your formula for frequency modulation implements a behavior very different to that of most synthesizers, and will sound uncanny to a musician.

For example, if lfo_y oscillates between -1 and 1 ; lfo_mod_depth is 0.5 ; and if note_frequency is equal to 220 Hz, osc_frequency will sweep between 110 and 330 - that is to say between -1 octave and +1 fifth around the note. So the pitch modulation will appear to be centered in the hertz scale, but will not be centered in the perceptual musical scale.

The correct behaviour is to have something like:

osc_frequency = note_frequency * 2 ** (tuning / 1200.0 + lfo_y * lfo_mod_depth)

Then:

tuning is expressed in cents, a musically relevant unit (100
cents = 1 semitone).
your pitch modulation is "musically centered",
and lfo_mod_depth is expressed in octaves.

This is called "exponential FM", and this is the norm on synthesizers. On analog synthesizers, this is implemented by summing the LFO signal to the CV that hits the exponential converter of the VCO. On digital synthesizers, this is implemented by applying the modulations on a high-resolution internal representation of pitch which is still on a musical scale - before converting to a frequency or a phase increment.


Different values when calculating signal's energy by hand vs MATLAB

Since I don’t know rest of the given code, I’ve decided to calculate the energy other way round.
The MATLAB code for calculating the given function’s energy is given below:
syms n; % Symbol assignment for summation

Energy_1 = symsum(power((2.3704 * exp(0.287682 * n)), 2), n, -Inf, -3); % Energy for the left-hand side of the sample domain

Energy_2 = symsum(power((2.3704 * exp(-0.287682 * n)), 2), n, 3, Inf ); % Energy for the right-hand side of the sample domain

Energy = Energy_1 + Energy_2 % Overall energy

According to this code above, the resultant overall energy becomes $E = 4.5715$.
As the number of samples increases, the result gets accurate. Also, the coefficient in the exponential’s power affects the accuracy of the result too, but this is not the main case.


When was the time-frequency analysis invented?

According to the preface of Foundations of Time-Frequency Analysis, a rough timeline is as follows:

1930 - Early development of quantum mechanics by H. Weyl, E.Wigner, and J. von Neumann.
1946 - Theoretical foundation of information theory and signal analysis by D. Gabor (cf. "Theory of communication").
1980 - Time-frequency analysis established as an independent mathematical field (apart from engineering) by Guido Janssen.
1990 - Development of wavelet theory. Overview of the mutual influence given in "The wavelet transform, time-frequency localization and signal analysis" by Ingrid Daubechies.


Can we broadcast two messages using NOMA?

There is a THRESHOLD, both can receive. If that would not have been possible why would NOMA be so popular, What is the use of NOMA then? Well, theoretically is possible however practically such channel state information is not available, that is what we read, however, I haven't seen whether that practically possible or not, Energy-Efficient NOMA Systems tells exactly what you have just said.


Why is each window/frame overlapping?

Windowing is a means to stationarize signals. Inside a small enough window, you can expect that the  properties of the signal chunk do not vary too fast. And you now can use tools well-suited to stationary signals, like Fourier-based techniques.
You can imagine non-overlapping rectangular windows, each defining a frame. Each sample is, somehow, treated with the same weight (the height of the window). However, when you look at the features extracted from two consecutive frames, or when processing them, the change of property between the frames may induce a discontinuity, or a jump ("the difference of parameter values  of neighbouring frames can be higher"). The same phenomenon happens with JPEG images, using $8\times 8$ non-overlapping blocks. This tiling or blocking effect is already annoying for the eye, it is even more disturbing for the ear.

Adding to this, if the signal evolves slowly with time, the left-end or the right-end sample could be seen as the "most different" samples: the frame is inherent localized at a specific location, depending on the window shape. If the window is symmetric, the central location is the center of the frame (some use asymmetric windows for speech processing).

So you can imagine giving some smaller weights to the edges of the frame. This is what you do when you are averaging a signal with unimodal weights, like triangular or Gaussian, instead of uniform. This can be beneficial for frequency estimation, as other windows have better behavior than the uniform one. However, this creates some imbalance in the processing (even information loss as said by @hotpaw2, esp. if the window vanishes) of some samples (the edge ones), unless you compensate it in some way, for instance by choosing a sum of windows whose amplitude ("constant overlap-add") or energy is close to one. Some figures can be found in
Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new at-top windows, and the framework of oversampled filter banks provides a quite unified treatment of such constraints.

My question is how the overlap lower/increase the processor demand?

It depends a lot of the implementation.  When the processing is done offline, you can save a lot of redundant operations, by folding the windows into independent blocks, subsampling where you can, etc. Using polyphase matrices or lifting from the overlapped filter bank theory can help you a lot in saving operations. Also, CPU and GPU aware processing can completely change the computational burden deal.
So what follows is quite crude.

Without much optimization however, if you have an overlap of $p$%, each sample will be processed, in average, $c = 1/(1-p)$ times. This is about the factor for the  number of frames that could be naively processed: $c=2$ for $p=50\%$, $c=4$ for $p=75\%$. If $p= 0$, frames can easily be processed in parallel. Otherwise, some memory or processing is required for the overlap.

At the extreme, you can perform a windowing centered at each sample. This is the principle at the basis of the time-frequency analysis. The computational cost can be high. But the high redundancy has a good side: any processing in the frequency domain would be more robust to noise, or imperfections in the processing: you should hear less clicks and pops.


Filter - Spatial Domain Versus Frequency Domain

It appears you need to study a bit on convolutional filtering of images, specifically on overlap add/ overlap save methods.

From the links I can see your objective is to apply the filter defined in the second equation i.e, the frequency domain one. The first link's accepted answer seems to be wrong (i didn't go through the code, clearly since the images shown by the paper and answer do not match).

Now the method you are using to apply the filter in the spatial domain is wrong. I'm assuming by padding you mean zero padding. The objective of zero padding before applying fft is to increase the resolution in the frequency domain. So imagine I have a 1D sine wave of 50 Hz. Taking fft of this signal would give a peak at the 50 Hz point with some additional values at different frequencies. This is because the sine is not pure sine since the signal is limited in time. If you had an infinitely long signal you would get a single value at 50 Hz. Padding the signal and then taking fft would still give the similar result but with more resolution on what additional frequencies exists.

In contrast convolution filter is essentially convolving 2 signals.
$$
y(t) = \int_{\tau = 0}^{\inf} h(\tau)x(t-\tau)d\tau
$$
In 2D, if we have a 3x3 filter kernel, we first multiply the first 3x3 block of the input with the kernel and then shift the kernel by one column. Then we multiply the next 3x3 with kernel and add the overlapping pixel values.

for small kernel sizes (say 3x3 or 5x5) this approach maybe faster. Using fft is alluring when working with large kernel sizes. The approach is go block by block and use the overlap save/add method (you can find a lot of info on this just google it for 2D). The reason we need to use overlap method is because fft facilitates circular convolution and not linear convolution (but circular convolution can be made to behave like linear convolution by adjusting the input and kernel sizes).

Once you understand the method then it is just deciding on the kernel size and using the equation for those u,v values. Hope this helps.


What's the noise in this signal ? (Beginner question)

From the spectrogram (frequency domain plot) you have a large signal at 60 Hz and harmonics. This will be mains pickup if this was recorded in a area with a 60 Hz mains supply. On the upper time plot this signal has a peak to peak amplitude of about 0.5 units. (Note that there are 30 cycles over 500 ms.)

It looks like to me the isolated peaks are the noise.

As JRE and MBaz suggest, it helps to know what your wanted signal looks like. The Wikipedia article on ECG gives an idealized view of a ECG trace (as well as some actual plots). You should be able to see that the isolated peaks are actually part of your wanted signal. Your patient seems to have a pulse rate of about 75 bpm.
All your wanted signal seems to be below 25 Hz. To get rid of the mains interference one way would be to low pass filter the signal or use a comb notch filter to remove 60 Hz and its harmonics. A notch filter may have less effect on your wanted signal but whatever filter is used you need to consider how its amplitude and phase response will affect your wanted signal.
An alternative approach is subtraction. You generate a local copy of the interfering signal and subtract it from the raw input leaving behind the wanted ECG. As this isn't based on a filter there are no unwanted amplitude or phase effects on the signal.


If you can use the Bilateral Filter then you can use the Guided Filter.
The nice property of the Guided Filter is its low complexity.
There is a simple and efficient implementation with with linear complexity of the number of pixels.
You may also have a look at:

Fast Anisotropic Curvature Preserving Smoothing.
Simple Image Edge Preserving Filter.


Any of the Stimson airborne radar books are good modern resources for general radar signal processing.
For more specific applications (basic STAP, EW, mismatched filters etc), you can look at Mark Richard’s book(s).
If your company has a library, check to see what radar books they have on hand; you may be able to get these for free.
Last, MATLAB has entire packages based around radar systems and plenty of examples that use them. If you have access to a MATLAB license, those might be good to check out as well.


The answer is clearly 'no', except for trivial cases. Note that
$$G(f)=|G(f)|e^{j\phi(f)}\tag{1}$$
where $\phi(f)$ is the phase of $G(f)$. Rewriting $(1)$ as
$$|G(f)|=G(f)e^{-j\phi(f)}\tag{2}$$
and taking the inverse Fourier transform gives
$$\tilde{g}(t)=(g\star a)(t)\tag{3}$$
where $\tilde{g}(t)$ is the inverse Fourier transform of $|G(f)|$, and $a(t)$ is the impulse response of a (generally non-causal) all-pass with frequency response $A(f)=e^{-j\phi(f)}$.
If we are given $|G(f)|$, we do not know $\phi(f)$ and, consequently, we don't know $a(t)$, so we can't compute $\tilde{g}(t)$ from $(3)$. Anyway, $(3)$ shows the relation between the inverse Fourier transform of $|G(f)|$ and $g(t)$.
Note that since $|G(f)|$ is real-valued, $\tilde{g}(t)$ must be even
$$\tilde{g}(t)=\tilde{g}(-t)\tag{4}$$


Compensating for group delay in FIR with symmetrical coefficients

No.

There is no way you can compensate for a  broad-band group delay without violating causality. If you add zeros at the end of your signal, the group delay stays constant and if you add zeros in front of it, the group delay goes up.

At some frequencies it may look that there is no group delay (or phase shift), but that's just an illusion: the phase of the filter at this frequency just happens to be an integer multiple of two pi. If you look at any other frequency, it won't be true any more.


Where does the "Desired Response" come from in a Wiener filter?

In many applications a Wiener filter is used to decorrelate the desired signal from the input to the filter. The desired system output is then the error signal, which is decorrelated from the filter's input signal.

As an example, in adaptive noise cancellation we assume that we have some noise reference signal (e.g., noise picked up by a microphone). And we have a noisy signal which we want to clean up. The noise in the signal is assumed to be correlated with the noise reference. In terms of Wiener filtering, the noisy signal is the desired signal, the noise reference is the input to the filter, the filter output is hopefully a good estimate of the noise in the desired signal, and the error signal is the actual output of the noise canceler, which is a cleaner version of the originally noisy signal. The figure below shows the set-up (from S.J.Orfanidis: Optimum Signal Processing):

There are many other applications of (adaptive) Wiener filters, such as adaptive equalizers in communication systems (see this answer), and echo cancelation (see this answer).

In sum, the "desired signal" is often a signal contaminated with noise or interference, and the output of the Wiener filter is an estimate of that noise or interference, which is subtracted from the "desired signal". This is possible if at the input of the filter there is some reference of the contaminating signal. The "error signal" is then the actually desired output, which is the noisy "desired signal" minus the noise-estimate at the output of the filter.


You have to estimate a velocity based on a regularly sampled signal. The velocity is given by the derivative on a continuous time signal. The approximation of this value with discrete data is a long-lasting topic. @Juancho is correct with the proposal of the $2$-point difference scheme.

As the sampling is never perfect, and $x[i]-x[i-1]$ is not symmetric. So if you do not need fully causal results, there exists more symmetric derivative schemes, such as an average of forward and backward differences:  $\frac{1}{2}(x[i]-x[i-1])+(x[i+1]-x[i]) = \frac{x[i+1]-x[i-1]}{2}$. Other schemes can be found in Numerical differentiation formulae (see here for a discussion).

If you need to compute forward derivatives for real-time purposes, a paper (in French, but formulae are readable) is Estimation par maximum de vraisemblance de la dérivée d'un signal bruité, with cares about noise issues.


Just sample a Sinc function at the appropriate phase offset, and truncate the FIR as needed to meet your rectangular window length requirement.  You only need 1 phase for a constant delay while keeping the same sample rate.


I don't understand how the zeros are found of this 2nd order feed-forward filter

Below is a polynomial in $z$ :
$$ P(z) = 1+\alpha_1z^{-1}+\alpha_2z^{-2} $$
Zeros of this polynomial is also known as its roots; i.e., those $z_r$ for which $P(z_r)=0$.
Solving for the roots of a general $n$-th order polynomial is not possible. They are numerically found.
However, one alternative way to define a polynomial is to write it as the product of factors of the form $1 - d_k z^{-1}$, where each $d_k$ is a zero/root of $P(z)$; i.e., $(1 - d_k d_k^{-1}) = 1-1 = 0$.
Then your second order polynomial $P(z)$ can be written as
$P(z) = (1-Z_1z^{-1})(1-Z_2z^{-1})$
where the complex valued $Z_1=Re^{j\theta}$ and $Z_2=Re^{-j\theta}$ are the roots of the polynomial. You should solve the second order equation $P(z)=0$ to find out what those zeros are.


I am not a specialist in image processing but I highly doubt this is feasible with your requirements (single image, no prior information). The challenges is twofold:

Detect which area of the image is a reflection. This looks like a quite complex scene understanding problem. In your example, I know that the right side is a reflection because I can infer from the image that a person is depicted here, and I have prior knowledge about what the back of a person looks like.
Reconstruct the background in the identified region. While it would be possible to just inpaint the area detected as a reflection with textures from the surrounding, the result would hide many details! Whatever I've seen in the field of single input image separation was relying on very simple image models (stationary textures), not real world images. The more traditional image separation setups (multiple channels, recto/verso of translucent pages, multiple polarization) do not apply to your case.


There are several problems in your code. First, it looks like you're confusing iteration and vector indices. The computation of e should use all values of the current (delayed) data vector, filtered with the current filter coefficients. In the update equation, you subtract a scalar from a vector, which is not what you want. Again you should be using all values of the current data vector (the length of which must equal the chosen filter length).

Take a look at the Matlab code in this question (the second one in the EDIT-part). I haven't run it but it looks like it deals correctly with the vectors in the update loop.


Frequency-selective channel: Basic in theory

In theory you need to recompute your precoding for every subcarrier.

In practice you might not want to do this as the computational complexity is pretty high. What one then typically does is grouping adjacent subcarriers and time indices together into "chunks" (I believe they are called resource blocks in LTE?) and just compute one precoding strategy for the entire chunk. The size of the chunk is a trade-off: if we make it too large, the quality of the precoding degrades because the channels becomes less and less similar over frequency and time. If we make it too small, computational complexity becomes an issue.

To determine the effective channel for the entire chunk that is used to compute the precoding, various strategies exist. The simplest might be just averaging the channels, for the MIMO case we often prefer to average the second-order statistics.


Question about delayed sampled sinusoid math expression

Indeed a typo
$e^{j \omega ((n - k)T)}$ would be the delay for $k$ samples, the $M$ in the expression $e^{j \omega(nT - M)}$ is the delay (in seconds if $\omega$ is in $rad/s$)


Minimum SNR Requirements for Maximum Likelihood Frequency Estimation

I doubt that there is such a rule because your assumption- the negative of computational work more than offsets the benefits of lower error probabilities- assumes that computations are expensive and errors are not. The "costs" of computational work and errors depends entirely on the application.


What does it mean for the DFT phase to be relative to a cosine wave?

A cosine wave is an even function around the origin (X axis = 0).  The basis vectors of the real components of a DFT are all cosine waves that start at 1 at the first sample.  Imaginary components of DFT basis vectors are all sine waves that are odd functions around the origin 0, and start at zero at the first sample in a DFT window.

However, for strictly real input, you can add a phase offset p to a cosine wave function to make cosine(x+p) look like any ratio mix of real and imaginary basis vectors (cosine waves and sine waves) in the DFT result, and that phase offset is what this book is calling the DFT phase.

For even length FFTs of non-fictional data, I prefer to do an fftshift to move the phase reference to the middle of the FFT window, instead of leaving the phase reference at the potential circular discontinuity at the first sample of the DFT window.


To simplify, let's consider the transmission of a single pulse $$s(t) = A p(t),$$ where $p(t)$ has energy 1, and the transmitted energy is $A^2$.

Note that the value of $A$ is related to the transmitted energy. If the transmitter's amplifiers have a lot of gain, then $A$ will be large. The AWGN channel model does not make explicit mention of this, but you should always keep it in the back of your mind that it's only a model, in a real system there is a lot more going on.

Now, the receiver has a filter matched to $p(t)$; when sampled at the right time, the matched filter output is equal to $A$. The receiver's noise has PSD equal to $N_0/2$, which is constant and independent of $A$.

You can then easily conclude that increasing the transmitted power will increase the SNR at the matched filter's output, which will in turn decrease the probability of error.

So, one possible approach would be:

Define a fixed transmitted constellation. Let's say you use BPSK with constellation $\lbrace -1, 1 \rbrace$.
Define a transmitter gain G.
Calculate the SNR at the filter's output: $2G^2/N_0$.
Plot the BER as a function of this SNR.


I've tried a copule of methods on your data (low-pass filtering via DFT, SSA filtering), but I don't know what result will satisfy you.

You mentioned that the signal is sampled with a higher-then-enough frequency. Can you rise the  sampling frequency even higher? I'm talking about oversampling. This can help to get rid of some part (some kind) of the noise present.

You also mentioned PCA and that it's not appropriate for univariate time series. Right, PCA is almost useles with univariate time series, but there is an appropriate and closely related technique, SSA. You can find detailed description here: Singular Spectrum Analysis for Time Series (free sample is enough: Chapter 2: Basic SSA). And my own GNU Octave (hope Matlab compatible) implementation here: ssa-octave.m (there is also Scilab implementation, but you probably don't need it). Take a look at chapter 2.4 "Choice of Parameters in Basic SSA" for interpretation and tuning tips.

It's also not clear to me: the "5 seconds of stationary" pattern in isometric data sample should be filtred out or it is a useful signal too?


there are a lot of misconceptions in the question.

"The received signal is an arbitrary waveform and not [only] a sine wave of frequency $f$."

it's a collection of sine waves at different frequencies.  and "tuning" a radio or a filter or whatever is selecting frequency components in the neighborhood of $f$ while rejecting the others.


Sensitiviy and Resolution in Windowing Functions

I object to the Wiki author's use of the words "resolution" and "sensitivity." Those words are vague and the author does not define or explain them in any meaningful way. If I was forced to use the author's poor terminology I'd say, "A fine-resolution window is a window having a narrow main lobe. (And windows whose main lobe widths are narrow have high-level sidelobes.) A window with good sensitivity is a window with low-level sidelobes. (And windows with low-level sidelobes have wider main lobes.)

I do not understand the 2nd and 3rd paragraphs of your question.

A window with a narrow main lobe improves your ability to distinguish (recognize, discriminate), in the frequency domain, two closely spaced spectral sinusoids.  A window with low sidelobe levels improves your ability to detect, in the frequency domain, a low-level sine wave that is close in frequency to a high-level sinusoid.


Why is PSNR used for image quality metrics instead of SNR?

Consider 8-bit quantization. This gives a mean squared error of about 0.577 least significant bit (LSB). Quantizing two images of different average intensities gives two different SNR, because that metric is normalized against the average intensity of the input image, but identical PSNR, because it is normalized against the maximum pixel value of 255. When describing the quality of something like a compression method, it is better to use a metric that is (to an extent) independent of the source material, because it carries meaning even when communicated without providing the input image.


how is this method of proving shift invariance correct

The so called example system is nonlinear and recursive and those books generally exclude the treatment of nonlinear & recursive systems. Hence you may not find example solutions.

The method relies on the fact that the relationship between the output $y[n]$ and input $x[n]$ remains the same between $y[n-d]$ and $x[n-d]$.

As an absurd example, think of the following. Consider the equation of a circle at origin: $x^2 + y^2 = r^2$. Now, shift this circle to the point $(a,b)$ and re-write its equation: $(x-a)^2 + (y-b)^2 = r^2$. Finally, redefine the variables $x' = x-a$ and $y' = y-b$ and rewrite the equation: $(x')^2 + (y')^2 = r^2$.

It's seen that the circle is shift invariant (not in the sense of the signal processing context though)


Adding "realistic" noise to signals of different amplitudes

I would do it this way in Matlab:

stanDevNoise_i=sqrt((Ai^2)/(2*EbNo));
noise_i=stanDevNoise*(randn(length(xi),1)+1i*randn(length(xi),1));

The reason I put (Ai^2)/2 is because you want to put there your signal power.
The noise is just a random number with the same length as your signal (p, in your case).


Interpolation techniques for extremely limited hardware

For retro hardware, the low pass filter after the PWM is probably the inertia of the speaker cone and/or your ear's internals. So first make sure the ringing isn't beyond your requirements at the full PWM sample rate, or due to quantization (use some tiny 20kHz sample rate waveform loops for testing).

For upsampling, I might try a zero or first order hold followed by an IIR filter (1st order or a biquad low-pass) chosen so that it only has a tiny number of "1"'s in the binary coefficients, allowing shift-and-add recursive filtering with deterministic latency.  If you can find coefficients with a tiny enough number of "1"'s in total, a cascade of biquad IIR's is also possible.

On the data preparation side, you can pre-noise-filter your samples on a more powerful system, instead of just depending on the retro-CPU's quantization.  You an also simulate your retro system's output filter, and tune the noise filtering to move some of the noise outside the passband.


Your second option is the most common: choose a frequency of interest, set the corresponding coefficient in the Fourier domain then take the inverse transform.
You have to pay attention to the correct scaling of the non-zero coefficient and to the particular FFT coefficient arrangement (continuous frequency in the middle or in one corner of the Fourier matrix).


Converting frequency from $\textrm{Hz}$ to radians-per-sample

Problems like these are best attacked using some dimensional analysis:

$$f_{[\rm rad/samples]} = f_{[\rm cycles/sec]}\cdot \frac{\text{sec}}{\text{samples}}\cdot \frac{\text{rad}}{\text{cycle}}$$

$$f_{[\rm rad/samples]} = f_{[\rm cycles/sec]}\cdot \frac{2\pi}{f_s}$$

where $f_s$ is the sample rate in $\textrm{Hz}$.


What is the difference between a Fourier transform and a cosine transform?

The Discrete Fourier Transform (DFT) and Discrete Cosine Transform (DCT) perform similar functions: they both decompose a finite-length discrete-time vector into a sum of scaled-and-shifted basis functions. The difference between the two is the type of basis function used by each transform; the DFT uses a set of harmonically-related complex exponential functions, while the DCT uses only (real-valued) cosine functions.

The DFT is widely used for general spectral analysis applications that find their way into a range of fields. It is also used as a building block for techniques that take advantage of properties of signals' frequency-domain representation, such as the overlap-save and overlap-add fast convolution algorithms.

The DCT is frequently used in lossy data compression applications, such as the JPEG image format. The property of the DCT that makes it quite suitable for compression is its high degree of "spectral compaction;" at a qualitative level, a signal's DCT representation tends to have more of its energy concentrated in a small number of coefficients when compared to other transforms like the DFT. This is desirable for a compression algorithm; if you can approximately represent the original (time- or spatial-domain) signal using a relatively small set of DCT coefficients, then you can reduce your data storage requirement by only storing the DCT outputs that contain significant amounts of energy.


There's no problem letting the coefficients be a slowly varying function of time,

$$S(t)=a(t)\sin(t)+b(t),$$

for example a linear ramp
$$a(t)=(a_1-a_0)\frac{t-t_0}{t_1-t_0}+a_0$$ for $t\in[t_0,t_1]$ and similar for $b(t)$.


How do delay lines create so many effects?

A DSP / Audio delay line doesnt create any effect in itself, a delay is simply a timed delay in the signal. it has no sound difference, it repeats the sound once and mixes it over original or sends it to a different mix line.

A feedback delay repeats the sound multiple times at the delay interval frequency, it sounds like a digital echo, and can become an oscillator above 20hz. almost all fx delays use variable feedback controls to control the number of repetitions of a sound, and they are filtered out as they feed back therefore diminishing.

Modern delays can work like grain synth samplers and often are based on them which have integrated effects like pitch shifting and research Grain-Delay and Grain-Resampling for infor about that. Squarepusher is somneone who uses alot of grain delay sometimes...

So a delay in effect is barealy an effect, but it can become an osciallator, has complex variants, and can feed into other lines that have any kind of additional FX.

To get a firmer grasp of that program i foudn this vid was more informative:

https://www.youtube.com/watch?v=ySzjFU6J4Ck

it's a pretty simple and very well balanced grain-delay with good controls and some bit crushing types of FX, and probably some filters too.


Relationship between $y(t)$ and $Y(f)$ in LTI system

The Fourier transform is linear, so a linear combination of multiple signals in the time domain also results into a linear combination in the frequency domain. So the Fourier transform of,

$$
y(t) = a_1\, x_1(t) + a_2\, x_2(t),
$$

yields the following relation in the frequency domain,

$$
Y(f) = a_1\,X_1(f) + a_2\,X_2(f),
$$

with $G(f)=\mathcal{F}\left\{g(t)\right\}(f)$.

However if instead you mean convolution, instead of addition, then you can use the fact that convolution in the time domain translates to multiplication in the frequency domain (the opposite is also true). So your example then becomes,

$$
y(t) = x_1(t)\, \large{*} \, \normalsize{x_2(t)} \equiv \int_{-\infty}^{\infty} x_1(\tau)\, x_2(t-\tau)\,d\tau,
$$

$$
Y(f) = X_1(f)\,X_2(f).
$$


How do I implement a nonorthogonal quadratic spline wavelet into discrete wavelet transform?

I would think it is available in the Matlab codes from the Wavelab toolbox. It hosts codes for papers and books. Especially, there is a directory \Wavelab\Books\WaveTour\ for Stéphane Mallat's book: A Wavelet Tour of Signal Processing. You can check images chapter by chapter from
A Wavelet Tour/Images. Those images can be regenerared. For instance, from Chapter 5,  Wavelab codes in \Matlab_Wavelab\Books\WaveTour\WTCh05 say from instance in wt05fig04.m:

  disp('The dyadic wavelet is a quadratic spline with one vanishing
  moment.') disp('These graphs display the modulus square of its Fourier
  transform') disp('dilated by 2^j for 0

which may suit your needs.


Recursive Least Square Adaptive Linear Equalizer

To answer (1) the adaptive equalizer without a training sequence (blind equalization) can be used based on the decisions of the received sequence. This specifically is called a "decision directed equalizer". Of course it can not work in very low SNR conditions, where a training sequence would be required. A typical approach is to have the training sequence first to remove channel ISI and establish a higher SNR signal. From this point a decision directed approach can continue to track and maintain the equalized signal as the channel varies with time. (If I recall correctly, a decision directed approach can typically converge if the bit error rate is better than 1E-3 -- but this exact figure I am not confident in).

I have heard (but have no experience with) that the most commonly used adaptive algorithm for blind channel equalization is the Constant Modulus Algorithm. Perhaps someone more knowledgeable about this approach can add an additional response with more comments and details about CMA (or better if I cannot Google a satisfactory explanation I can post it as a question myself).

To answer (2), the LMS (Least Mean Square, also referred to as the Gradient Algorithm) and the RLS (Recursive Least Squares algorithm) are the two most common algorithms for recursively minimizing mean square error in adaptive linear equalizers (in contrast to non-linear Decision Feedback Equalizers, which are not to be confused with "Decision-Directed" mentioned above, would be preferred in the case of frequency selective fading).

The LMS algorithm is well known, easy to implement and computationally cheap (requires $2M+1$ multiply operations, where M is the number of equalizer coefficients). It uses a Stochastic Gradient Descent Rule, so may converge slowly.

The RLS algorithm in comparison is computationally more expensive (requires $2.5M^2+4.5M$ multiply operations), but has faster convergence and good tracking.

For details on the LMS and RLS algorithms (and likely where I got the metrics quoted above), see Wireless Communications by Rappaport which I have notated in the plots below. Notice the computational difference between the two algorithms:

Also see this post here where I give the code and detailed intuitive explanation for the determination of the LMS equalizer coefficients using the Wiener-Hopf equation: Compensating Loudspeaker frequency response in an audio signal. In this post I used a known copy of the signal to determine the channel impulse response in the static case (not a recursive solution). However as noted above this could also be done blind with the decisions of what we think was transmitted if the SNR was high enough to make a sufficient number of correct decisions.

you can find additional detailed explanations on Wikipedia: https://en.wikipedia.org/wiki/Recursive_least_squares_filter


Does DFT produces the same output as FFT?


Difference Fourier Shift Theorem and Cross Correlation Theorem

First of all, what you quote as the Cross correlation Theorem

  the cross-correlation between two signals is equal to the product of fourier transform of one signal multiplied by complex conjugate of fourier transform of another signal

doesn't make much sense: the product mentioned above gives the _cross_power spectral density function, and not the cross-correlation function at all. The cross-correlation function is the inverse Fourier transform of the cross-power spectral density function, namely, the product mentioned above.

That being said, the key difference between what you are doing and what the Cross-Correlation Theorem is that you are using a special case that is perhaps better termed the Autocorrelation function theorem.  Notice the use of the word another (next to last word in the highlighted wording above). When another is replaced by the same, we get that the autocorrelation function is the inverse Fourier transform of the power spectral density.


Even though the signals are sampled you can get accuracy which is well above the accuracy offered by the samples as long as you sample using Nyquist.
Actually, Using the Matched Filter you can achieve the CRLB (Cramer Rao Lower Bound) for Delay Estimation (Easy to derive for white noise).
If you calculate the CRLB for Time Delay Estimation you'll see it depends on the Signal BW and the SNR and not the sampling frequency (Given Nyquist sampling).
How can that be achieved?
You can interpolate the the cross correlation signal and by simple math infer the time delay in resolutions well beyond the sampling resolution.
For instance, with signals with BW of ~ 15 MHz we shall sample by 60 Mhz.
The sampling resolution (Assuming RADAR) is ~ 2.5 [Meter] yet the CRLB at decent SNR will give you a bound of few Centi Meters which can be achieved using the Matched Filter + Interpolation.
Enjoy.


Haar Training with object Blacked out in Negative Samples

I have had good success doing this with my training sets. I have modified the OpenCV_Annotation application to present a fixed rectangle that moves with the mouse, once the space bar is pressed the coordinates are written to the positives.txt file and a copy of the image is saved into a pre-defined negatives subdirectory after graying out the rectangle.
I have seen no training problems or performance issues although I use other negatives as well as the blurred out images. The application uses a gray background to overwrite the identified area instead of just blacking it out. I am unsure if that would make a difference in the negative.


Zero-frequency response of an IIR transfer function

Because of the pole near the unit circle the impulse response doesn't really settle until after tens of thousands of samples. In comparison, your input of 15 samples is very short. Try making it longer:

x = rand(15E6,1);
y = filter(b,a,x);
mean(y)/mean(x)

>> 0.34252

Also, the filter starts in a state that corresponds to an all-zeros input history. If you precondition it with a run of the actual input, the result will converge faster:

x = rand(15E6,1);
y = filter(b,a,x);
mean(y(5E6:length(y)))/mean(x(5E6:length(x)))

>> 0.34483


A feature extraction method: Principal Component Analysis (PCA)

Ok, it sounds like you are trying to do eigenfaces, right?

In that case, you have to think of your face images as points in a very high-dimensional space. For example, if your images are 32x32, then the space has 32 * 32 = 1024 dimensions.  Operating in so many dimensions is very difficult, because distances between points become almost meaningless.  With so many dimensions any point is very far away from any other point, even if the two points represent the same face.

That's why you want to reduce the dimensionality, and PCA is one way to do it. If you go from 1024 dimensions to, say, 10, then the distances will become more meaningful, which fits in with your idea of "highlight their similarities and differences".  However, PCA by itself will not tell you of any patterns in your data, beyond the directions of greatest variance.  PCA will project your data into a lower-dimensional space with decorrelated axes, but you have to do more work to actually discover the patterns in the data, such as train a classifier, or run a clustering algorithm.

Now, to understand what happens to the edges and transitions from light to dark, you have to display your eigenvectors as images and look at them. You will see that you don't really preserve the precise edge information, but you do keep course-level information about light and dark areas.


amplitude of upsampled and downsampled signal without filter

Close, but not the right formula for the up-sampler:

$$x_i[n] = x[n/L]$$

The "up sampler" is really an expander that inserts zeros between the samples of x[n], and it has the formula:

$$
x_u[n] =\begin{cases}
x[n/L] & n=0, \pm L, \pm 2L, ... \\
\\
0 & otherwise
\end{cases}
$$

The DTFT of the up-sampler "expander" can be derived as follows:

using $\delta(n)$ we can rewrite $x_u[n]$ as:

$$
x_u[n] = \sum_{k=-\infty}^{\infty}x[k]\delta[n-kL]
$$

Applying definition of DTFT to $x_u[n]$:

$$
\begin{aligned}
X_u(\omega) &= \sum_{n=-\infty}^{\infty} \left(x_u[n] e^{-j\omega n} \right) \\
\\
X_u(\omega) &= \sum_{n=-\infty}^{\infty} \left( \sum_{k=-\infty}^{\infty}x[k]\delta[n-kL] \right) e^{-j\omega n} \\
\end{aligned}
$$

Applying sifting property of $\delta$:

$$
\begin{aligned}
X_u(\omega) &= \sum_{n=-\infty}^{\infty} x[k] e^{-j\omega L k} \\
\\
X_u(\omega) &= X(\omega L)
\end{aligned}
$$

Thus, the gain of up-sampler "expander" in freq domain is unity.


How to apply the provided wave transforms on an image

You want to apply to your original image $I_1[n,m]$ a transform (x,y)=T(n,m) such as a warping effect, to obtain a new image $I_2[n,m]$.

In your example, you want to shift original pixel positions (n,m) to new coordinates based on the calculated position (x,y) from the transform equation (x,y) = T(n,m). As usually the case, this direct transform of indices will possibly result in noninteger positions (x,y) and therefore be quantized to slightly distorted locations, or even worse some pixels of the final image may not be filled, leaving them empty.

One approach to solve these two problems is to apply the algorithm from the reverse, by using the inverse transform of T(,)  provided it exists. As in the pseudo matlab code here: Result is as shown below

As a specific case of your examples:
WAVE1: T(n,m) = (n + 20*sin(2*pi*m/128) , m ) then the inverse transform is
       invT(x,y) =( x - 20*sin(2*pi*y/128) , y) which yields m = y and n =  x - 20*sin(2*pi*y/128)

% NOTE1: This code is not efficient. You shall Vectorize it easily
% and get rid of the for-loops, yielding a substantial boost in efficiency.
% NOTE2: boundary rows/columns may be left untouched which needs a simple fix
% The transform that you want to appply is (x,y)=T(n,m)
% The transform that we will use is its inverse (n,m)=invT(x,y) (you shall invert your transform)
% Original image I1 is of size S1xS2, transformed image I2 is of size N1xN2
for x=1:N1    % let N1 = S1 and N2= S2 for this example WAVE1
  for y=1:N2

     % s1- compute back coordinates in the original image domain
     % [n m] = invT(x,y);
     m = y;
     n = x - 20*sin(2*pi*y/128);

     % s2-s3 Apply (bilinear) interpolation to find the image intensity, if
     % n,m does not fall into integer indices
     d1 = n - floor(n);
     d2 = m - floor(m);

     % s3- define the image pixel
    if ( ((n>=1)&(n<S(1))) & ((m>=1)&(m<S(2))) )
     I2(x,y) = (1-d1)*(1-d2) * I1(floor(n),floor(m)) + ...
               (1-d1)*d2 * I1(floor(n),floor(m)+1) + ...
               d1*(1-d2) * I1(floor(n)+1,floor(m)) + ...
               d1*d2 * I1(floor(n)+1,floor(m)+1) ;
    end%if

  end
end


I'm interested in SSA and the related techniques and I made an implementation which may be useful to you. Take a look at this Scilab code: https://gist.github.com/werediver/9785544 Scilab language is pretty similar to Matlab and the code should be easy to port.

Now to your questions and your code.

You're applying SVD directly to the trajectory matrix (the time delay embedding of the time series). That's possible and will work, but usualy SVD is applied to the covariance matrix $C=XX'$.

Right after SVD you have the eigenvalues (they are shown on your plot), but SSA is not done yet. Next, you need to compute the principal components of the trajectory matrix: $V_i=\frac{X'U_i}{\sqrt{\lambda_i}} (i=1,...,L)$ You need $\sqrt{\lambda_i}$ when applying SVD to the covariance matrix and just $\lambda_i$ when SVD is applied directly to the trajectory matrix.

So now the SSA decomposition is done and you want a partial reconstruction. The time domain representation of $V_i$ is $X_i=\sqrt{\lambda_i}U_iV_i'$. These components are additive so sum all the components you're interested in and then apply diagonal averaging to get from $X_I$ to univariate time series. Diagonal averaging produces a vector from averaged antidiagonals of the $X_I$ matrix, take a look at my implementation for example.

I can point you to a good description of Basic SSA: Singular Spectrum Analysis for Time Series book by N. Golyandina & A. Zhigljavsky has freely available sample, Chapter 2: Basic SSA. This chapter completely covers Basic (univariate) SSA decomposition & reconstruction algorithms and gives many advices on interpretation of the decomposition.

There are also pretty good book (teaching aid) in Russian language if it's appropriate for you: Метод «Гусеница»-SSA: анализ временных рядов  by N. Golyandina.


Based on your information a packet length is about 77.2 millisecond. With 3918 bps, each packet contains about 302 bits. Frequency is shifting about every 0.255 millisecond. you get about 256 samples per bit, so 128 point FFT seems to be correct. I assume you have a real signal, so you need to evaluate only 64 bins. If you are synchronized properly, each two consecutive FFTs will provide you with one frequency, switching at the other. If you are not synchronized with the start of the first bit, only every other bin will provide you with a nice single frequency (assuming that you do not have noise and spurious signals)


1D Deconvolution with Gaussian Kernel (MATLAB)

Approaches

There are many methods for Deconvolution (Namely the degradation operator is linear and Time / Space Invariant) out there.
All of them try to deal with the fact the problem is Ill Poised in many cases.

Better methods are those which add some regularization to the model of the data to be restored.
It can be statistical models (Priors) or any knowledge.
For images, a good model is piece wise smooth or sparsity of the gradients.

But for the sake of the answer a simple parametric approach will be takes - -Minimizing the Least Squares Error between the restored data in the model to the measurements.

Model

The least squares model is simple.
The objective function as a function of the data is given by:

$$ f \left( x \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

The optimization problem is given by:

$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

Where $ x $ is the data to be restored, $ h $ is the Blurring Kernel (Gaussian in this case) and $ y $ is the set of given measurements.
The model assumes the measurements are given only for the valid part of the convolution. Namely if $ x \in \mathbb{R}^{n} $ and $ h \in \mathbb{R}^{k} $ then $ y \in \mathbb{R}^{m} $ where $ m = n - k + 1 $.

This is a linear operation in finite space hence can be written using a Matrix Form:

$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| H x - y \right\|_{2}^{2} $$

Where $ H \in \mathbb{R}^{m \times n} $ is the convolution matrix.

Solution

The Least Squares solution is given by:

$$ \hat{x} = \left( {H}^{T} H \right)^{-1} {H}^{T} y $$

As can be seen it requires a matrix inversion.
The ability to solve this adequately depends on the condition number of the operator $ {H}^{T} H $ which obeys $ \operatorname{cond} \left( H \right) = \sqrt{\operatorname{cond} \left( {H}^{T} H \right)} $.

Condition Number Analysis

What's behind this condition number?
One could answer it using Linear Algebra.
But a more intuitive, in my opinion, approach would be thinking of it in the Frequency Domain.

Basically the degradation operator attenuates energy of, generally, high frequency.
Now, since in frequency this is basically an element wise multiplication, one would say the easy way to invert it is element wise division by the inverse filter.
Well, it is what's done above.
The problem arises with cases the filter attenuates the energy practically into zero.  Then we have real problems...
This is basically what's the Condition Number tells us, how hard some frequencies were attenuated relative to others.

Above one could see the Condition Number (Using [dB] units) as a function of the Gaussian Filter STD parameter.
As expected, the higher the STD the worse the condition number as higher STD means stronger LPF (Values going down at the end are numerical issues).

Numerical Solution

Ensemble of Gaussian Blur Kernel was created.

The parameters are $ n = 300 $, $ k = 31 $ and $ m = 270 $.
The data is random and no noise were added.

In MATLAB the Linear System was solved using pinv() which uses SVD based Pseudo Inverse and the \ operator.

As one can see, using the SVD the solution is much less sensitive as expected.

Why is there an error?
Looking at a solution (For the highest STD):

As one could see the signal is restored very well except for the start and the end.
This is due to the use of Valid Convolution which tells us little on those samples.

Noise

If we added noise, things would look differently!
The reason results were good before is due to the fact MATLAB could handle the DR of the data and solve the equations even though they had large condition number.

But large condition number means the inverse filter amplify strongly (To reverse the strong attenuation) some frequencies.
When those contain noise it means the noise will be amplified and the restoration will be bad.

As one could see above, now the reconstruction won't work.

Summary

If one knows the Degradation Operator exactly and the SNR is very good, simple deconvolution methods will work.
The main issue of deconvolution is how hard the Degradation Operator attenuates frequencies.
The more it attenuates the more SNR is needed in order to restore (This is basically the idea behind Wiener Filter).
Frequencies which were set to zero can not be restored!

In practice, in order to have stable results on should add some priors.

The code is available at my StackExchange Signal Processing Q2969 GitHub Repository.

This answer was given to a different question yet both of them deals with deconvolution of a signal with a Gaussian Kernel.

If there are issues not answered, please point me and I will broaden.


How can I classifying this License Plate images? Matlab , LPR ,

Here's a procedure you should try in a first attempt :

You have RGB images, so each images is a 3*RC matrix, where R is the number of rows and C the number of columns.
We will consider each pixel of an image as a realization of a random vector $\mathbf{x}$ of dimention 3. Hence, having an image, you have RC realization of $\mathbf{x}$. With the classification exemple you gave, each classe has 4 images, so we have 4*R*C realization of $\mathbf{x}$. We are abble to estimate some stastical moments, that we hope will be good descriptor, such as mean, covariance etc ... (mean can be estimated with empirical mean or median).
Let's just stick with the mean for now. Let us call $m_1,m_2,m_3,m_4$ the means of each class you gave.
Now we have a new image we want to classify, we estimate it's 3-dimensional mean $\mathbf{\mu}$, then using nearest neighbour, the class is given by :
$$
C=\arg_i\min\left(\lVert\mathbf{\mu}-m_i \rVert _2\right)
$$

This is very basic.
You may improve this scheme using additional features like variance (diagonal of the covariance matrix) or covariance matrix. You may estimate variance or covariance using empirical variance or empirical covariance estimator (realy just the classical formulae). However the euclidian distance won't be a good distance since you will comparing inhomogenous vector (the features will be composed of mean + covariance), you'll want to use a more advanced distance such as Bhattacharyya distance for the special case of gaussian distribution.