The role of GPS in INS/GPS navigation systems

By itself, no, it cannot.

Your method of fixing the INS system is pretty clever, but it only accounts for drift in the gyroscope. Accelerometers are also (very) prone to drift, generally showing up as an increasing velocity bias. The typical method to account for accelerometer drift is to use zero-velocity updates. By detecting points when you are not moving, you can periodically remove any non-zero velocities from your model.


Phone conversation diarization with LIUM

system.Telephone is just for detection of speech, that's why you get only init clusters.  You can use default fr.lium.spkDiarization.system.Diarization it will do proper steps as described in LIUM docs, result will be like this:

;; cluster S0
file 1 0 1135 M S U S0
file 1 5722 300 M S U S0
file 1 9266 359 M S U S0
file 1 12670 370 M S U S0
;; cluster S1
file 1 1135 832 F S U S1
file 1 5372 350 F S U S1
file 1 6022 572 F S U S1
file 1 10167 1550 F S U S1
file 1 13040 466 F S U S1
file 1 13623 352 F S U S1
;; cluster S10 [ merge HCLR 0 = S10 in S11 with 1.208107727685253 ]
file 1 1967 693 F S U S10
file 1 2660 1644 F S U S10
file 1 4304 1068 F S U S10
file 1 6594 937 F S U S10
;; cluster S12
file 1 7531 429 F S U S12
file 1 8692 574 F S U S12
;; cluster S13
file 1 7960 732 M T U S13
;; cluster S16
file 1 9625 542 F S U S16
file 1 11717 953 F S U S16

LIUM models included into jar are trained for 16khz, so the input should be better 16khz, you need to resample your 48khz file to 16khz before you feed into LIUM.

If you want to run on telephone speech, you might need to retrain mail/female/ubm models for that. It is also worth to note that LIUM models are trained for French, for other languages they should be better retrained.

Overall, LIUM is not expected to work out of box, you might need to spend some time to make it work properly.


Why are equalizers moving average model (FIR)

The reason why almost all linear adaptive equalizers are implemented as FIR filters is that FIR filters are always stable and that there exist relatively simple and effective adaptation algorithms. Note that much work has been done on adaptive IIR filters (e.g., this book by Phillip Regalia), but in practice FIR filters are still the preferred option.

Note that no real-world channel can be ideally equalized by an FIR filter, but that's not as bad as it may seem because an FIR filter can approximate the ideal equalizer reasonably well if its filter length is chosen appropriately and if the adaptation algorithm is tuned to the application. Also, perfect equalization of the channel is not always desired if the channel is not only dispersive (i.e., acts a filter) but if it also introduces noise. In this case a combined criterion is used that takes into account the ISI due to the dispersive channel as well as the degradation due to the noise.


Do Nyquist samples have the same sum as the integral of the signal?

The equation in your question is correct. If we assume that $x(t)$ is an ideally band-limited signal with bandwidth $B$, then, according to the sampling theorem, we have

$$x(t)=\sum_{n=-\infty}^{\infty}x(nT)\frac{\sin\left[\pi (t-nT)/T\right]}{\pi (t-nT)/T}\tag{1}$$

where the sampling rate $1/T$ satisfies $1/T>2B$.

If we assume that the integral over $x(t)$ exists, using $(1)$ we obtain

$$\begin{align}\int_{-\infty}^{\infty}x(t)dt&=\int_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}x(nT)\frac{\sin\left[\pi (t-nT)/T\right]}{\pi (t-nT)/T}dt\\&=\sum_{n=-\infty}^{\infty}x(nT)\int_{-\infty}^{\infty}\frac{\sin\left[\pi (t-nT)/T\right]}{\pi (t-nT)/T}dt\\&=\int_{-\infty}^{\infty}\frac{\sin(\pi t/T)}{\pi t/T}dt\cdot \sum_{n=-\infty}^{\infty}x(nT)\\&=T\sum_{n=-\infty}^{\infty}x(nT)\tag{2}\end{align}$$

where the last equality follows from the fact that the integral in the one but last line just equals the DC value of an ideal low pass filter with gain $T$.


How can the I and Q components of a quatrature signal be different from each other (apart from the 90º phase shift)

You are mistaken in your 2nd paragraph.  The I and Q components are not always identical with a 90 degree phase offset.  Instead, they are usually generated by multiplying some real or complex signal with two sinusoids that offset by 90 degrees.  Only if the input is a constant (e.g. carries zero information) does the IQ output fit your description.

A strictly real signal has no imaginary component.  But when heterodyne modulated (complex multiplied) by a complex carrier (two higher frequency sinusoids, cosine and sine), the result is a complex IQ signal.  If the baseband signal is not just strictly real, but complex, the result of that heterodyne modulation can produce completely independant and unrelated I and Q components in the IQ result.

When down-converting a higher frequency strictly real signal to baseband, again you will get identical I and Q with a 90 degree offset only for constant input.  But if the real input is not constant, but varies over time, the complex multiplication result can produce very different I and Q component waveforms.  This is because the local maxima of the cosine waveform and sine waveform occur at different points in time, and thus can be multiplied by different input values if the input is varying, thus producing in different result values, thus making the waveforms different, not just offset by 90 degrees.

Since QAM can carry information (e.g. not just a single constant value), the I and Q components can vary.


You can easily verify yourself that the following two ways of cascading two filters are equivalent (up to numerical inaccuracies):

% method 1
y1 = filter(h1,1,x);
y = filter(h2,1,y1);

% method 2
h = conv(h1,h2);
y = filter(h,1,x);

So you were on the right track when you convolved the two impulse responses. Note that linear convolution is what you need, not circular/cyclic convolution.

As for the multiplication with $z^{-q}$, where I hope that $q$ is an integer, you need to understand that it is not the (time-domain) signal that is being multiplied, but that this multiplication refers to the Z-transform of the signal. In the Z-transform domain, a multiplication with $z^{-q}$ corresponds to a delay of $q$ samples.

Finally, regarding your expectation that the cut-off frequency of the cascaded system is $\pi/4$, note that the frequency response of the cascaded system is simply the multiplication of the two individual frequency responses, so the cut-off frequency does not change significantly (it may change a bit depending on how you define it). The transition from pass band to stop band will be steeper, and the stop band attenuation increases, but the start of the transition band remains unchanged. If both filters are identical you will basically get a better low pass filter with the same pass band and stop band edges.


What is the main difference between the acoustic and electromagnetic waves?

Reciprocity, particularly in Air acoustics.  If you are down wind from from an acoustic source like an animal and you hear it, you are less likely to be heard by it.  The same effect is true with currents in underwater acoustics but to a lesser extent because sound velocity is relatively greater than fluid velocity.

The variation of sound speed tends to be more important in simulations. There are specializations of the wave equation, actually more typically the Helmholtz equation to account for things like $c(z)$, depth dependence.

Acoustic waves do not have polarization, but do have vector particle velocity.

At the other extreme, the speed of light can not be exceeded while the speed of sound isn’t a universal constant.

The Friis equation is not used in acoustics.  Isotropic lossless propagation in acoustics is usually inverse R squared. The term “path loss” have different conventions. The Friis equation implicitly includes the aperture of the receiver, so for isotropic lossless propagation, the path loss is proportional to frequency.  Radio engineers often (stubbornly) claim that low frequencies propagate further because of the Friis equation. Lower frequencies do often propagate further but typically it is because refraction is related to wave length, so in terrestrial point to point communications, the curve of the earth, scatter by objects like houses, penetration into buildings, and mode coupling between the earth’s surface and ionosphere favor low frequency propagation.  High frequencies tend toward line of sight.  Friis isn’t why.

In acoustics, performance is typically ambient noise limited while in radio,  the limit is usually related to the antenna temperature. When an array is moving through water, flow noise can exceed ambient noise. The calculation of array gain differs accordingly.

Acoustic waves are often encountered in combination with sheer waves.  There are also slow waves in EM propagation but sheer wave propagation is essential in seismology.

In biological systems, first principle bulk properties like the Young’s modulus of living lung tissue tend to be avoided because they are very complex to model and directly measure, particularly if you don’t want to kill the subject.

There is also a distinct area of research relating acoustics to perception.  Humans (and animals) are better at localizing sounds than classic array theory predicts.


How Last Edge Image Can Be Achieved from Law Masks

If I understand correctly, the question is, given many images which are result of different Edge Filter applied on the same image, how to actually mark edges.

Well, you basically created 25 tests for each pixel to decide whether or not it is an edge.

You could apply many approaches to decide:

Majority Votes - If more than half of the voters decided it is an edge, nark it as an edge.
Threshold - If more than $ x $ voters vote for edge the pixel will be declared as an edge.
Spatial Model - Instead of per pixel decision, look around it and other voters.
Weight of Votes - Don't mark votes as "Yes / No" but give it a scalar. If the sum of all scalars above a threshold value, declare an edge.

As you can see the options are endless.


Filters generally can't just have a single cut-off frequency (otherwise their complexity will be huge).  So a transition region is usually needed.

In this case, the filter_stop_freq is that frequency below which the filter MUST act like a stop filter and filter_pass_freq is that frequency above which the filter MUST act like a pass filter.

The frequencies between filter_stop_freq and filter_pass_freq are the transition region or band.

Because this is designing an FIR filter, filter_order is just the length (number of taps) of the filter. The higher the order, usually, the better the fit (and the smaller the transition region can be).


Conditions on $v \in \mathbb{C}^{n}$ such that the complex to real DFT is sensible

Obviously any operation is allowed that keeps the symmetry of the frequency domain data intact, i.e. any operation that doesn't add an imaginary part to the time domain sequence. And this is the case for all standard filtering operations unless you want to generate a complex valued (time-domain) signal by some operation. This could be the case if you want to generate an analytic signal, but in that case you actually make two signals out of one (the real part and the imaginary part in the time domain). So for normal low pass, band pass, high pass etc. filtering you don't need to worry. As long as you do things right, the symmetry of the time domain sequence should not be changed by these operations.


Researchers from the Johns Hopkins University have recently released a corpus of music, speech, and noise which, according to them, is suitable for training models for voice activity detection and music/speech discrimination.

See https://arxiv.org/pdf/1510.08484.pdf for details.


Estimate Instantaneous Frequency Using LMS Algorithm

I make no claim this is optimal and I also do not think that the overshooting you mention is necessarily incorrect. Nothing can change instantaneously, so you are bound to see some range of frequencies when you transition from one sinusoid to the next.

The first idea I had for a least-squares estimate was to use an itegrator on the unwrapped phase. Since frequency is the derivative of the phase, if we integrate the frequency we will get the phase. In terms of a linear system, this looks something like

$$
\mathbf{Ax} = \boldsymbol{\phi},
$$

where

$$
\mathbf{A} =
\begin{bmatrix}
1  & 0 & 0 & \cdots \\
1 & 1 & 0 & \cdots \\
1 & 1 & 1 & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{bmatrix}
$$

and $\boldsymbol{\phi}$ is the unwrapped phase. To get the unwrapped phase, you need to make the real-valued signal analytic by putting it through an I/Q demodualation process. Then the least-squares solution for $\mathbf{x}$ in the above will be an estimate of the frequency.

The octave code I cobbled together is below. Since the problem is small, I explicitly form $\mathbf{A}$, but for larger problems you would want to use a function handle with something like lsqr() in order to solve the system. I also run a moving average over the frequencies so that the jump from one to another is a bit smoother, which should give less overshooting.

randn( 'seed', 71924 );
rand( 'seed', 7192 );

fmin = 100;
fmax = 1e3;
fs = 1.2 * 2 * fmax;
T = 5;
t = ( 0 : 1/fs : T ).';
Ns = length(t);

% Define some noise.
snr = 10.^( 50 / 10 );
n = sqrt( 1/snr ) * randn( Ns, 1 );

% Choose some random frequecies.
Nf = 50;
f = fmin + ( fmax - fmin ) * rand(Nf,1);

% Choose the intervals.
tmp = randperm( Ns );
inds = sort( tmp( 1 : Nf - 1 ) );

% Define the frequency for each sample.
fvec = zeros( Ns, 1 );
last = 1;
for ii = 1 : length(inds)
  fvec( last : inds(ii) ) = f(ii);
  last = inds(ii) + 1;
end
fvec(last+1:end) = f(end);

% Moving average over the frequencies to prevent
% bandwidth when the frequency changes.
fvec = conv( fvec, 1/3 * ones(3,1), 'same' );

% Define the sinusoid with time-varying frequency.
% We do so in a way that makes it phase continuous.
% We also add the noise here.
x = cos( 2*pi * cumsum( fvec ) ./ fs ) + n;

% The frequency sample locations.
df = fs / Ns;
floc = ( 0 : df : ( fs - df ) ) - ( fs - mod( Ns, 2 ) * df ) / 2;

% Shift the center of the positive half of the
% spectrum to zero.
y = x .* exp( -1j * 2*pi * t * fs/4 );

% Design a half-band filter.
Ntap = 51;
N = Ntap - 1;
p = ( -N/2 : N/2 ).';
s = sin( p * pi/2 ) ./ ( p * pi + eps );
s( N/2 + 1 )= 1/2;
win = kaiser( Ntap, 6 );
h = s .* win;

% Form an analytic signal.
z = conv( y, h, 'same' ) .* exp( 1j * 2*pi * t * fs/4 );

% The unwrapped angle.
phi = unwrap( angle( z ) );

% Use an intergrator to get a least-squares
% estimate of the phase.
A = tril( ones( Ns, Ns ) );
fest = ( A \ phi ) * fs / (2*pi);

figure();
hold on;
plot( t, fest );
plot( t, fvec );
legend( 'Estimated', 'Truth' );
xlabel( 'Time (s)' );
ylabel( 'Frequency (Hz)' );


How to find PDF value from histogram plot in MATLAB?

The histogram is a good numerical estimate for the PDF sampled a the bin centers. Just scale it so that the sum over all bins is 1.
The quality of the estimate depends on number of samples and "good" choice of the bins. The latter one is a trade off: smaller bins give you better amplitude resolution and wider bins make the PDF in each band more accurate.
In your case this looks fairly uniformly distributed. The two outer bins are lower than the rest which could be real, but it could also be caused by the bins extending past the max and the min of a signal. It's important how to line up the bins with the min, mean & max of the signal.


High FIR order in time critical classification task

I would handle this as a chain of FIR filters with decimation to get much lower tap count FIR filters that don't need to be nearly as steep in the early stages.
Maybe using Blackman-Harris windowed filters for the first two stages (for superior alias fold-in rejection).
Then a Kaiser windowed filter for the 3 stage to get a flatter top and a sharper roll-off.
I guarantee you'll be using much less than 6000 taps total for the same filtering performance.
I've worked out fairly decent scheme in some Octave/MatLab code:
LPF -> decimate by 2 -> LPF -> decimate by 5 -> BPF
This gets you to a final 100 samples/second sample rate.
The three filters, that I came up with, have tap lengths of 9, 51, and 205 respectively.
The two low-pass anti-alias filters keep any alias fold-in into your band of interest at -40 dB or lower.
The final band-pass filter is -35 dB at DC, -6.3 dB at 0.5 Hz, -0.1 dB at 22.5 Hz, and -40 dB at 24.5 Hz.

Here's the Octave code.
Main script:
% First anti-alias filter, designed for decimation by 2 after the filter
Fs1 = 1000;  % [-500,500]
a = blackman_harris_lpf(1.0, Fs1, 100.0, 300.0, 55);
length(a)

% Second anti-alias filter, designed for decimation by 5 after the filter
Fs2 = Fs1/2;  % [-250, 250]
b = blackman_harris_lpf(1.0, Fs2, 50, 20.0, 45);
length(b)

% Band pass filter for the frequency band of interest
Fs3 = Fs2/5;  % [-50, 50]
c = kaiser_lpf2(1.0, Fs3, (23.5-0.5)/2, 2, 9);
radian_phase_inc = pi/(Fs3/2) * (24.5-0.5)/2;
c = 2*real(c .* exp(1j * radian_phase_inc .* [0:(length(c)-1)]'));
length(c)

% Get response curves for each filter for plotting
[ha,wa] = freqz(a, [1], 2*5*1000);
[hb,wb] = freqz(b, [1], 2*1000);
[hc,wc] = freqz(c, [1], 2*1000/5);

% Create plotting x-axis (aka frequency-axis) value arrays to
% fold freq axis to show aliasing into final passband.
% Scale to original Fs1
waturn = [wa(1:(2*5*1000/2))' wa((2*5*1000/2):-1:1)'] * 500/pi; % Fs = 1000, fold for decimation by 2
wbturn = [wb(1:(2*1000/5))'   wb((2*1000/5):-1:1)' ...
          wb(1:(2*1000/5))'   wb((2*1000/5):-1:1)' ...
          wb(1:(2*1000/5))'                         ] * 250/pi; % Fs = 500, fold for decimation by 5
wcturn = wc * 50/pi; % Fs = 100, no turns

clf;
plot( ...
     waturn, 20*log10(abs(ha)), ';1st filter, Blackman-harris LPF, 1000 sps folded to 500 sps;', ...
     wbturn, 20*log10(abs(hb)), ';2nd filter, Blackamn-harris LPF,  500 sps folded to 100 sps;', ...
     wcturn, 20*log10(abs(hc)), ';3rd filter, Kaiser BPF, 100 sps;'  ...
);
title('Filters Alias Folding Performance');
ylabel('Gain (dB)');
xlabel('Frequency (Hz)');
grid on;
hold off;

blackman_harris_lpf.m:
function b = blackman_harris_lpf(gain, Fs, Fc, W, atten_dB)
% B = blackman_harris_lpf(GAIN, FS, FC, W, ATTEN_DB)
%
% Design a Blackman-harris windowed low pass filter
%
%    GAIN Passband linear gain. E.g. 1.0
%
%    FS   Sampling frequency in Hz
%
%    FC   Cutoff frequency (the edge of the ideal LPF response before
%         before truncation and windowing of the sinc() pulse) in Hz.
%
%    W    Transition width in Hz
%
%    ATTEN_DB  stop band attenuation in dB

M = floor(atten_dB * Fs/(22.0*W)); % rule of thumb from fred harris
if mod(M,2) == 0  % ensure M is odd
   M = M + 1;
end

% Windowed, truncated sinc()
w = blackmanharris(M);
h = sinc([-(M-1)/2:1:(M-1)/2]*Fc/(Fs/2)).' .* w;

% Normalize passband response to 0 dB and apply gain
b = gain * h / sum(h);

kaiser_lpf2.m:
function b = kaiser_lpf2(gain, Fs, Fc, W, Beta)
% B = kaiser_lpf2(GAIN, FS, FT, W, BETA)
%
% Design a Kaiser windowed low pass filter
%
%    GAIN Passband linear gain. E.g. 1.0
%
%    FS   Sampling frequency in Hz
%
%    Fc   Cutoff frequency (the edge of the ideal LPF response before
%         before truncation and windowing of the sinc() pulse) in Hz.
%
%    W    Transition width in Hz
%
%    BETA Kaiser window parameter dictating the width of the main lobe
%         and the height of the stopband lobes.  Beta == 0 is equivalent
%         to a rectangular window.

M = floor((Beta/0.1102 + 8.7) * Fs/(22.0*W)); % rule of thumb from fred harris
if mod(M,2) == 0  % ensure M is odd
   M = M + 1;
end

% Windowed, truncated sinc()
w = kaiser(M, Beta);
h = sinc([-(M-1)/2:1:(M-1)/2]*Fc/(Fs/2)).' .* w;

% Normalize passband response to 0 dB and apply gain
b = gain * h / sum(h);


Most of the classical linear transformations (even filtering) may have three main types of scalings:

natural or none: sum or integral does not have an explicit scaling factor (but the inverse may need one)
in amplitude: because of linearity, at least one reference "unit" signal should have a unit amplitude in the transformed domain, for an easy read of graphs. Typically a unit amplitude sine (which is not $L_1$-integrable) should have a one amplitude at its frequency in the Fourier domain
in energy: to have constant energy in the original or the transformed domain.

Each of these options (sometimes, two of them are equal) depends on the purpose. Fourier  has a specific relation to scale (Where in “Discrete Time Signal Processing” (Oppenheim et al.) can I find the scale-change theorem?):
$$  s(\alpha t) \mapsto  \frac{1}{|\alpha|} S \left(  \frac{f}{|\alpha|}\right) $$
Obtaining scale-, shift- or orientation-invariant features is a long-lasting issue in DSP, and is still a concern in machine intelligence or deep learing.
Along the Fourier line, you can look at the Fourier-Mellin transform or scale represensation (by L. Cohen), and the many feature representations like SIFT, SURF, ORB, BRIEF: Image Matching Using SIFT, SURF, BRIEF and
ORB: Performance Comparison for Distorted


Take the Z Transform of your difference equation and perform algebraic manipulations to find $H(z)=\dfrac{Y(z)}{X(z)}$ as a rational function of $z$.

Take the Z Transform of your input $x(nT) = x[n]$ to find $X(z)$.

Multiply to get $H(z)X(z) = Y(z)$.

Find the inverse Z Transform of $Y(z)$ to find $y[n]$.  (This will be the difficult part, as you may need to use partial fraction expansion and other mathematical techniques so you can look up inverse transform in a table.)

These tables will help:
http://lpsa.swarthmore.edu/LaplaceZTable/LaplaceZFuncTable.html

Update in response to comment

Here's what I get.

Finding the transfer function:

$$\mathcal{Z}\left\{y[n]\right\} = \mathcal{Z}\left\{y[n-1] + x[n] -x[n-K]-x[n-L]+x[n-K-L]\right\} $$

$$Y(z) = Y(z)z^{-1} + X(z) - X(z)z^{-K} -X(z)z^{-L} + X(z)z^{-K-L}$$

$$\begin{align*}H(z) = \dfrac{Y(z)}{X(z)} &= \dfrac{1-z^{-K}-z^{-L}+z^{-K-L}}{1-z^{-1}}\\
\\
&= \dfrac{\left(1-z^{-K}\right)\left(1-z^{-L}\right)}{1-z^{-1}}\\
\\
&= \left(1 + z^{-1} + \dots + z^{-K+1}\right)\left(1-z^{-L}\right)\\
\\
&= 1 + z^{-1} + \dots + z^{-K+1} -z^{-L}-z^{-L-1}- \dots -z^{-K-L+1}\\
\\
\end{align*}$$

Transforming the input:

$$X(z) = \mathcal{Z}\left\{x(nT)\right\} = \mathcal{Z}\left\{Ae^{-\frac{nT}{\tau}}u[n]\right\} = A\dfrac{z}{z-e^{-\frac{T}{\tau}}}= A\dfrac{1}{1-e^{-\frac{T}{\tau}}z^{-1}}$$

Filtering in the Z domain:
$$\begin{align*}Y(z) &= H(z)X(z)\\
\\
&=\dfrac{1-z^{-K}-z^{-L}+z^{-K-L}}{1-z^{-1}} \cdot A\dfrac{1}{1-e^{-\frac{T}{\tau}}z^{-1}}\\
\\
&=A\left[\dfrac{1-z^{-K}-z^{-L}+z^{-K-L}}{\left(1-z^{-1}\right)\left(1-e^{-\frac{T}{\tau}}z^{-1}\right)}\right]\\
\\
&=A\left[\dfrac{C}{1-z^{-1}}+\dfrac{B}{1-e^{-\frac{T}{\tau}}z^{-1}}\right]\\
\\
\end{align*}$$

Performing partial fraction expansion:

$$\begin{align*}1-z^{-K}-z^{-L}+z^{-K-L} &= C\left(1-e^{-\frac{T}{\tau}}z^{-1}\right) + B\left(1-z^{-1}\right)\\
\\
\left(1 + z^{-1} + \dots + z^{-K+1}\right)\left(1-z^{-L}\right)\left(1-z^{-1}\right)&=C\left(1-e^{-\frac{T}{\tau}}z^{-1}\right) + B\left(1-z^{-1}\right)\\
\end{align*}$$

so
$$\begin{align*}C &= 0 \\
\\
B &= 1 + z^{-1} + \dots + z^{-K+1} -z^{-L}-z^{-L-1}- \dots -z^{-K-L+1}\\
\end{align*}$$

Continuing filtering in the Z domain:
$$\begin{align*}Y(z) &=A\dfrac{1 + z^{-1} + \dots + z^{-K+1} -z^{-L}-z^{-L-1}- \dots -z^{-K-L+1}}{1-e^{-\frac{T}{\tau}}z^{-1}}\\
\\
\end{align*}$$

Since there is no information on the values of $K$ and $L$, there is not really any more simplification we can do in the Z domain.

Taking the inverse Z-Tranform:

$$\begin{align*}y[n] &= A\left[e^{-\frac{nT}{\tau}}u[n]+e^{-\frac{(n-1)T}{\tau}}u[n-1]+\dots+e^{-\frac{(n-K+1)T}{\tau}}u[n-K+1]-e^{-\frac{(n-L)T}{\tau}}u[n-L]-e^{-\frac{(n-L-1)T}{\tau}}u[n-L-1]-\dots-e^{-\frac{(n-K-L+1)T}{\tau}}u[n-K-L+1] \right]\\
\\
&= A \sum_{k = 0}^{K-1} e^{-\frac{(n-k)T}{\tau}}u[n-k] - e^{-\frac{(n-k-L)T}{\tau}}u[n-k-L]\\
\\
\end{align*}$$


Is There a Way to Perform a 2D Image Rotation by Matrix Multiplication?

Let's say you have an $ M \times N $ image.
If you turn it into a vector in $ {R}^{\left( M N \right)} $ and create a matrix which is $ \left(M  N \right) \times \left( M  N \right) $ by multiplication and rearranging the vector from the multiplication result you can reposition any pixels where ever you want.

Being practical, by padding the original image you can create rotation by any angle you want even with interpolation (By summing few pixels into one position).

It all starting with vectorization of the image.


You are asking 2 questions:

How to extract those labeled frequencies and remove the rest of the data?
How to do Demodulation?

Well, regarding your first question, a dirty solution would be by working the the DFT data.
Just zero all frequencies but those you're after and apply Inverse DFT (ifft).
Leave some margin around your wanted frequencies.

About the second question, look here:
http://en.wikipedia.org/wiki/Demodulation

Basically, in order to shift data in the frequency domain you multiply by pure harmonic signal in the time domain.
You chose the frequency of this harmonic signal s.t. (At least this is the classic case) the frequency you want will move to the DC frequency.
Then you apply Low Pass Filter to leave only the data you wanted.

This is the opposite of the Modulation process.


What you have for your 2nd 'resonance' is actually a vibrational mode with a transfer function of the form ('s' being the Laplace variable)

H(s) = (s^2 + 2*zeta1*omega1*s + omega1^2) / (s^2 + 2*zeta2*omega2*s + omega2^2)

The denominator 'omega2' term corresponds to the peak at 4.99GHz, and the numerator 'omega1' term corresponds to the deep trough that follows the peak at 4.99 GHz.

One can demonstrate this behavior in Matlab very simply using the code below

zeta1 = .01; omega1 = 1.02; zeta2 = .011; omega2 = 1.0;

num = [1 2*zeta1*omega1 omega1^2];

den = [1 2*zeta2*omega2 omega2^2];

freqs(num,den)

The frequency response is shown in the figure below which mimics your 2nd 'resonance'.

What you need to do for your 2nd 'resonance' is to estimate 4 (not 2) parameters, viz., omega2 and zeta2 (for the peak) and omega1 and zeta1 (for the trough). Depending on your data a nonlinear least squares optimization scheme should do. As a start, you could assume zeta1 and zeta2 to be equal (say zeta) and get an initial estimate of the terms 'omega1' and 'omega2' from your graph (or data). This leaves you with solving for just 'zeta'. Then you could iterate the model to get zeta1 and zeta2, etc., etc.

Note that the 'Q' term that you used is directly related to the quantity 'zeta' in my exposition here.


Preprocess sound for voice recognition in Pocket Sphinx on Raspberry Pi 3

Human voice has the same frequencies as other noises, it is not possible to filter it with dumb algorithm. It is possible to filter based on temporal patterns like spectral subtraction, but this way you can only filter static noise, not background speech/bable noise. Pocketsphinx already has spectral subtraction.

To filter noise reliably you need an advanced noise filtering neural network or something like that, see RNNNoise, it is very resource-intensive. Also, it is better to recognize noise sound instead of filtering it since recognizer can figure out things by himself.

One important approach is to reliably filter is to use direction of arrival of the sound, that is a microphone array which can filter noises based on direction. An example of such array is Respeaker or Matrix One. Thats why all smart speakers have microphone arrays as a critical component.

For more accurate recognition it is better to use modern toolkit with neural network like Kaldi, it is also much more accurate with noise filtering and without it.

Further reading for you:

Snips - Snips is a partially proprietary framework, but you can get some nice ideas and hardware advices. It is easy to get a good system using their packages too. Also see Community projects
Make a smart speaker with Respeaker
Zamia Speech the easiest way to run Kaldi on Raspberry Pi


Symbol table: $1+1j, 1-1j, -1+1j, -1-1j$

The Symbol table needs to be a python list/array - [1+1j, 1-1j, -1+1j, -1-1j]


Is a Kalman filter suitable to filter projected points positions, given Euler angles of the capturing device?

Low Pass Filtering

It would be good to know what you mean by "simple low pass filter".

For example, if your measurements at time $k$ are

$$
p_{k} = \left[ \begin{array}{c} x_k \\ y_k \end{array} \right]
$$

and your low pass filtered estimates are:

$$
p^{\tt LPF}_k = \alpha p^{\tt LPF}_{k-1} + (1-\alpha)p_k
$$

then you will have quite a large group delay in the filter of about $1/(1-\alpha)$ (for alpha close to 1).

Modeling the Signal: Simplistic Approach

To use the Kalman filter (or any similar approach), you need to have a model for how your measurements are acquired and updated.

Usually this looks like:

$$
p^{\tt TRUE}_{k+1} = {\mathbf A} p^{\tt TRUE}_k + {\mathbf B} \epsilon_k
$$
where $\epsilon_k$ is the process (driving) noise, ${\mathbf A}$ is the state transition matrix, and ${\mathbf B}$ is your input matrix.

And then your measured $p_k$ are:
$$
p_k = {\mathbf C} p^{\tt TRUE}_k + {\mathbf D} \nu_k
$$
where $\nu_k$ is the output (measurement) noise,  ${\mathbf C}$ is the output matrix, and ${\mathbf D}$ is your measurement noise matrix.

Here, the "state" of the model is chosen as the true positions, and the things you measure are the output.

You can then apply the Kalman filter equations to this to get state estimates $\hat{p^{\tt TRUE}_k }$ of the true position.

However, this approach is simplistic because it does not use any knowledge of how the points might move (nor does it use your 4 points and any knowledge you might have about how they move together).

Modeling the Signal: Starting a Better Approach

This page shows how to set up the problem involving the positions and euler angles.
It's doing something different from what you need, but the state is:

$$
p^{\tt TRUE}_{k} = \left[
x_k\  y_k\  z_k\
\dot{x}_k\  \dot{y}_k\  \dot{z}_k\
\ddot{x}_k\  \ddot{y}_k\  \ddot{z}_k\
\phi\ \psi\ \theta\
\dot{\phi}\ \dot{\psi}\ \dot{\theta}\
\ddot{\phi}\ \ddot{\psi}\ \ddot{\theta}\
 \right]^T
$$

and the measurements (output) is

$$
p_k = \left[ x_k\  y_k\  z_k\  \phi\ \psi\ \theta\ \right ]^T
$$

All the model on that page is really doing is saying:
$$
x^{\tt TRUE}_k = \sum_{n=0}^k \dot{x}^{\tt TRUE}_n n\Delta t + \frac{1}{2}  \sum_{n=0}^k \ddot{x}^{\tt TRUE}_n (n\Delta t)^2
$$
(but for each of $x,y,$ and $z$).

This is just the classic "equations of motion". See equation (3) here.


An "integrate and dump" filter is a common approach for this application, which is essentially an averaging filter averaged optimally over your symbol period. To implement this, you simply reset your accumulator at the symbol boundaries, accumulate your output over the symbol period and then "dump" your result by taking the accumulated output as the result before resetting it again. For a rectangular pulse shape BPSK in the presence of additive white Guassian noise, this is the optimum matched filter.

Note, unless you have a shared clock between your transmitter and receiver, you will likely also need to implement some form of carrier recovery to track clock variations between your transmitter and receiver (otherwise once you drift 90 degrees, your cosine will be a sine and you will get no output in that condition!). And as described above, you will need some form of timing recovery to establish when your start of the symbols are referenced to your signal as received.

I have detailed these in the following posts:

Phase synchronization in BPSK

Recovering signal for psk

Gardner Timing Recovery for Repeated Sybmols


A strong LOS component, i.e. a high Rician factor, will always improve receiver performance. If $K$ denotes the Rician factor, i.e. the ratio between the power of the LOS component and the power of the other (scattered) components, then for $K\rightarrow\infty$ you have no multipath propagation at all and you get a simple AWGN channel. For $K=0$ you have a Rayleigh fading channel (with no LOS component).


Noise is definetly undesired, but not incomprehensible.

The fact that noise is comprehensible is the reason we talk about it. Cosmic Rays might be present, but they have no effect on an old mechanical clock. Would you consider them "noise" when you read the clock?

You might call that loud disturbing sound that you hear "noise". A deaf person would not consider that noise, because it's incomprehensible.


How to sample continuous signal correctly?

close all;
clear all;

interval        = 0.5;
f0              = 16;
t               = 0:1/2048:interval-(1/2048)
x_t             = sin(2*pi*f0*t);

figure(1);
subplot(2, 1, 1);
plot(t, x_t);
axis([0 interval])
hold on;

fs              = 1024;
N               = (interval * fs);
sample_index    = [0:N-1];
x_n             = sin(2*pi*(f0/fs)*sample_index);

n               = 0:1/fs:interval-(1/fs)

stem(n, x_n);

X               = fftshift(fft(x_n, N));

df              = fs / N;
sampleIndex     = -N/2:N/2-1;
relative_f      = sampleIndex * df;

subplot(2, 1, 2);
stem(relative_f, abs(X));
xlabel('f (Hz)'); ylabel('|X(k)|');
axis([min(relative_f) max(relative_f)])


Is it possible to convert from 2D to 3D form?

[EDITED] Here's how it's done.

Steps:
1. Isolate the Road Divider Part.

Then, using Houghlines, find out the longest lines in Image. Find out the extrema points that cross image boundary. You got the Quadilateral points. I skipped this part by Manually Choosing them. In my case, the width of road at top of image is 10, and at bottom is 60.

Now, for the Destination, we have to Image so that Divider Quadilateral has uniform width of 60, so we save in another array new points, which have same bottom points, but top points are modified as they form a rectangle with height same as image (or bigger, you can adjust as it suits you after running it once and seeing the results), but width 60. Remember, we are transforming our image such that the divider fits into this Rectangular Area. Image will change automatically.

Now, using
Mat TransformMat = getPerspectiveTransform(ipPts, opPts);

of opencv on initial and destination quadilateral Points, we get a Transform Matrix, Which we apply to our Test Image.

warpPerspective(ipImg, opImg, TransformMat, ipImg.size());

You'll see a lot of Area is Cropped to fit the Output in InputImage Size. However you'll notice that Image is Transformed to get parallel Divider, as we Expected. By performing some Translation and Stuff, and taking a bigger Size Output Image, You'll get this. We can never get a perfectly wrapped Output Image, since its size can be very large. But we'll reach some constraint that'll fulfill our task. Here's Last image :


Frequency measurement from acoustic emission sensor data(VDC RMS)

Yes an FFT approach will work to identify frequency content, but you may need longer data sets to see something useful. For the FFT you will see the frequencies from 0 to your sampling rate (25000), and there will be the same number of samples in frequency as you have in time.  Therefore your frequency resolution is 25000/N where N is the number of samples.

Note, importantly the frequency resolution is dependent only on the time length of your data, T and is independent of the sampling rate. The above relationship works out as if you sample more, N will be proportionally larger in a fixed block of time T.

Therefore

$$ Res = 1/T $$

Where Res is the frequency resolution of your measurement and T is the time length of the data.

For your two datasets, you have 122 data samples for the first one, and and 222 samples for the second; which is 4.42 ms and 8.84 ms of data. The frequency resolution is 25000/122 = 204.92 Hz for the first block data and 25000/222 = 112.61 Hz for the second block (or 1/4.42 ms and 1/8.84 ms).

Let's look closer at your data, the second set that was larger. From a plot in the time domain (see below), we see a very large transient at the initial part of your capture that would likely obscure our ability to see the general frequency content, so I chose to work with a subset of your data from 2.8ms to 8.84 ms (which would decrease the frequency resolution further unfortunately but allows us to not be overcome by the initial impulse).

Showing the region after the initial impulse:

Here is the FFT of that portion of your sequence. I had removed the mean before taking the FFT which is why there appears to be a large null in the center.

The peak is near 500 Hz, with around 170 Hz of resolution.

This is what the FFT looks like without removing DC or the initial transient in the beginning.

Windowing your data first can help reduce the side-lobes in frequency to give you more dynamic range (if you are looking for lower level signals that may be masked by the side-lobes from larger signals), but know that that will decrease your frequency resolution further.

Here is the same plot as above of your entire second data set after applying a Kaiser window.

And a zoom in of above showing the FFT bin sample locations:

Zooming in further, we see the detail below. The Kaiser window used was with a Beta of 8, which has a frequency resolution that is approximately 3x wider than the non-windowed example. Since I used the complete data set (and the effects of the impulse are masking the lowest frequency region accordingly), the resolution non windowed would be 112.6 Hz, and with the window used around 330 Hz.

My recommendation to get more meaningful results is to increase the time length of your capture (assuming your frequency is stationary for a longer time). Given the frequency roll-off visible from these FFT's, it is not necessary to sample at 25 KHz; so you could reduce your sampling rate if that helps you increase the time length of the capture.


The ultimate book on the topic is probably http://www.amazon.com/Spatial-Hearing-Revised-Psychophysics-Localization/dp/0262024136 . It's a bit dated but contains and excellent review of the main psycho-acoustic research on spatial hearing.

It is a very complicated topic. For example azimuth and elevation perception work completely different and hence azimuth resolution is quite a bit better than elevation resolution. The actually resolution (or jnd for just noticeable difference) depends heavily on the actual and the properties of the signal.


No, scrambling does not increase the bandwidth.  Scrambling does even out the bandwidth, such that the power is spread evenly over the signal bandwidth.

You are thinking of chipping, which is used in DSSS systems.  Scrambling and chipping are different processes.


The problem is likely to be your filter. Designing a very narrow filter is very hard; a very large order is required and even then, some responses may simply be unrealizable. Note that "narrow" here is relative to the Nyquist frequency. I don't know what your sampling frequency is, but let's assume 8000 Hz. Run this code:

fn=4000;f1=59;f2=61;w1=f1/fn;w2=f2/fn;
b=butter(5,[w1 w2],'bandpass');
freqz(b);

The command freqz plots the actual frequency response of the filter that was returned by butter. It is very important always to verify that your filters' response resemble what you had in mind -- that is, that they meet your requirements.

I got slightly better results using fir1 instead of butter. However, in your case, what I would try to do is to downsample the signal so that, relative to the Nyquist frequency, the filter is not too narrow. To see the effect, run this code:

fn=120;
f1=59;f2=61;w1=f1/fn;w2=f2/fn;
b=fir1(200,[w1 w2],'bandpass');
freqz(b)

You'll see a nice, narrow response with 0 dB gain at the center frequency and around -70 dB rejection out of band.


This is a classic exercise and it shows you how to obtain the DFT of a real-valued length $N$ signal in terms of the DFT of a length $N/2$ complex-valued signal.

You have to write $X[k]$ in terms of the even and odd samples of $x[n]$. With $W_N=e^{-j2\pi /N}$ you get

$$\begin{align}X[k]=\sum_{n=0}^{N-1}x[n]W_N^{kn}&=\sum_{n=0}^{N/2-1}x[2n]W_{N/2}^{kn}+W_N^k\sum_{n=0}^{N/2-1}x[2n+1]W_{N/2}^{kn}\\&=\text{DFT}_{N/2}\{y_R[n]\}+W_N^k\;\text{DFT}_{N/2}\{y_I[n]\}\\&=Y_e[k]-jW_N^kY_o[k]\tag{1}\end{align}$$

where $y_R[n]$ and $y_I[n]$ are the real and imaginary parts of $y[n]$, respectively, and $Y_e[k]$ and $Y_o[k]$ are the even and odd parts of $Y[k]$. Eq. $(1)$ shows that after having computed the length $N/2$ DFT $Y[k]$, you need to split it into its even an odd parts, and from them you can obtain the length $N$ DFT $X[k]$.

Obtaining $Y[k]$ from $X[k]$ is also possible, even though that's usually the less practical problem, because often you have an FFT routine that takes complex input signals, so it's more efficient to compute a lenght $N$ FFT of a real-valued sequence by computing a length $N/2$ FFT of a complex-valued sequence. Anyway, $Y[k]$ is given by

$$Y[k]=\sum_{n=0}^{N/2-1}x[2n]W_{N/2}^{kn}+j\sum_{n=0}^{N/2-1}x[2n+1]W_{N/2}^{kn}\tag{2}$$

With $W_{N/2}^n=W_N^{2n}$, the first term on the right-hand side of $(2)$ can be written as

$$\begin{align}\sum_{n=0}^{N/2-1}x[2n]W_{N}^{k(2n)}&=\sum_{n=0}^{N-1}x[n]W_N^{kn}\frac{1+(-1)^n}{2}\\&=\frac12\sum_{n=0}^{N-1}x[n]W_N^{kn}+\frac12\sum_{n=0}^{N-1}x[n]W_N^{nN/2}W_N^{kn}\\&=\frac12\left(X[k]+X[k+N/2]\right)\tag{3}\end{align}$$

where I used $(-1)^n=W_N^{nN/2}$.

Eq. $(3)$ shows that decimation in the time domain corresponds to aliasing in the frequency domain, as expected.

In a completely analogous manner, the last term in Eq. $(2)$ can be written as

$$j\sum_{n=0}^{N/2-1}x[2n+1]W_{N/2}^{kn}=\frac{j}{2}\left(X[k]-X[k+N/2]\right)W_N^{-k}\tag{4}$$

Substituting $(3)$ and $(4)$ into $(2)$ gives the final solution:

$$Y[k]=\frac12X[k]\left(1+jW_N^{-k}\right)+\frac12X[k+N/2]\left(1-jW_N^{-k}\right)\tag{5}$$


Discrete time aperiodic signal: Concept of frequency?

TL;DR

Your statements are inaccurate. The concept of frequency is $\frac{1}{T}$ where $T$ is the period time. Be careful with the continuous-time and discrete-time definition for a perioud.

  discrete time sinusoids have frequency −π<ω<π or $−\frac{1}{2}<f<\frac{1}{2}$

So, this statement is not accurate.

A frequency is defined as the number of times a signal repeats itself in a second and is measured in $[Hz]=[\frac{1}{sec}]$. This is, by definition, not limited to the interval you stated and, infect, not defined well for negative values at all.

An angular frequency, $\omega=2\pi f$ is the rate of change of angular displacement, θ, (during rotation), or the rate of change of the phase of a sinusoidal waveform. It is measured in angle per second or $\frac{radians}{sec}$. As can be seen from the units, this is also not limited to an interval.

When we start to analyze using the Fourier and Laplace transforms, we analyze in the frequency domain. First, the Laplace transform analyzes the signal and if we want to consider frequencies we convert it to continuous-time Fourier transform (CTFT) by using $H^F(\omega)=H^L(s=j\omega)$, which exist iff $0\in ROC$ for the Laplace transform. Otherwise, CTFT is not defined.

CTFT is a function of $\omega$ and is not limited to an interval. However, let us look at a $T$ periodical signal $x(t)=x(t+T)$ with a CTFT $X^F(\omega)$.
$$x(t)=x(t+T)$$
$$\mathcal{F}\{x(t)\}=\mathcal{F}\{x(t+T)\}$$
$$X^F(\omega)=X^F(\omega)e^{-j\omega T}$$
$$X^F(\omega)-X^F(\omega)e^{-j\omega T}=0$$
$$X^F(\omega)(1-e^{-j\omega T})=0$$

Either $X^F(\omega)=0$, or $(1-e^{-j\omega T})=0$ and sinc $(1-e^{-j\omega T})=0\implies\omega T=2\pi k$ we know that $X^F(\omega\neq\frac{2\pi k}{T})=0$. This means that the CTFT of a periodic signal is discrete.

  Is the concept of frequency to a-periodic signal whether it is continuous or discrete-time is valid?

You have to be careful here. Is the discrete signal $\sin[n]$ a periodic signal? The surprising answer here is no. A descrit-time periodic signal is defined as a signal for which exist $N\in \mathbb{N}$ S.T. $x[n]=x[n+N]\forall n$. Let us consider $0=\sin[0]=\sin[N]\implies N=2\pi k \notin \mathbb{N}$. on the other hand $\sin(t)$ is most definitely periodical. The concept of frequency in both cases is $\frac{1}{T}$ where $T$ is the time of a single period. For the continuous-time, it is the smallest $T$ for which $x(t)=x(t+T)$ and for a periodic discrete signal, the smallest $N$ for which $x[n]=x[n+N]$. The frequency of both is $\frac{1}{T},\frac{1}{N}$ respectively.

the discrete-time Fourier transform (DTFT) is the analogy of CTFT for sampled\descrete signals. The difference between both is that it is periodical with respect to the sampling rate. The time distance between every two samples is $T_s$ hance the sampling frequency is $F_s=\frac{1}{T_s}$ and $\omega_s=\frac{2\pi}{F_s}=2\pi T_s$. we now take the original Fourier transform and replicate it around $k\omega_s$ where $k\in\mathbb{N}$. After this operation, we normalize (divide) by $\omega_s$. Now, the DTFT is a function of $\theta$ which is limited to $-\pi \leq \theta\leq \pi$ and is periodic around $2\pi$.
All this can be written as: $$X^f(\theta)=\frac{1}{T_s}\sum_{k=\infty}^{\infty}X^F(\frac{\theta-2\pi k}{T_s})$$

Maybe this is what you meant when you wrote the limitation on $\omega$?

I am not sure where did you get the final relation. Maybe a reference will help and I will edit the answer.


Is the system $y(n)=x(n^2)$ time invariant?

The system is obviously time-varying. However you make the following usual mistake in the first case of showing the response to the delayed input $x_2[n] = x[n-k]$ as follows.

Writing the output $y[n]$ in the formal manner: $$ y[n] = \mathcal{T}\{x[n]\} = x[n^2] $$ and denoting the shifted input as $$x_2[n] = x[n-k]$$ you would have the corresponding output $y_k[n]$ as:

$$y_k[n] = \mathcal{T}\{ x_2[n]\} = x_2[n^2] = x[n^2 -k] $$ which is defintely not equal to $y[n-k] = x[(n-k)^2]$ hence the system is time-varying.

You should have defined the output as $$ x_2[n] = x[n-k] \implies x_2[n^2] = x[n^2 -k] $$ but you didn't. See that the equality sign between $x[n]$ and $x_2[n]$ enables you to operate on $n$ only, not the whole argument of both sides.


Detecting defects/differences in printed images

All the comments in the comment section were very constructive, thank you.

In the end, I ended up using the solution here

For my case, most of my images are not expected to tilt more than a few degree, lighting and depths are strictly controlled. So a Euclidean Motion model (MOTION_EUCLIDEAN) should suffice for findTransformECC which computed the transformation matrix that is required to correct the misalignment.


Reconstruction of a signal using an ideal DAC

This is an ill posed question.

First it doesn't talk about the output DAC reconstruction sampling rate which is necessary for the reconstruction filter design and which also determines the frequency spectrum scaling of the reconstructed signal. Anyway then you may assume the DAC sampling-rate to be equal to ADC sampling rate for simplicity, which are both $F_s = 100$ Hz in your case. (btw your matlab code is wrong here, you have taken the sampling period to be $0.001$ which should be $0.01$.)

Furthermore, it talks about an ideal DAC which cannot be realized by any code, but, its output can be mathematically described perfectly without any computation by definition of the ideal DAC to be $$y_r(t)=x_c(t)$$ as long as there is no aliasing during the ADC stage which is the case here. So it's not clear if they want you to actually compute an approximate output $y_r(t)$ or simply evaluate the exact described output $x_c(t)$ in terms of input samples, which is a null computation.


What does "frequency" mean for various kinds of signals?

The key insight that Fourier had when he developed Fourier analysis is that any absolutely integrable (thanks Jason R) function can be represented as the weighted sum of sines and cosines.  Explaining why this is true is way beyond the scope of this answer.  I suggest you study Fourier theory to understand this better.


Calculate the entropy of Turkish alphabet and compare with the entropy of English alphabet with using Matlab?

I'll give you a hint: You need the probability distribution of the alphabet (English and Turkish). For English see the Wikipedia "Letter Frequency" page. Now it's possible to compute the Entropy of this distribution.

Note that this is a rough approximation (and an upper bound) since in practice, the letters of a word (and a sentence) are correlated. So, for example, in English the letter 'q' is always followed by 'u'. Thus, the entropy of a source producing English documents is strictly less than the one you'll get from looking at single letters.


Amplitude Response at greater than half the sampling frequency

This is a phenomenon called aliasing

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

Usually, there are anti-aliasing filters that will filter out signals with frequencies higher than fs/2.

However, in your case I think there isno anti-aliasing filter. Hence your 14 kHz-signal aliases to 2 kHz. Since f > fs/2, then apparent frequency = 16 kHz - 14 kHz = 2 kHz


In the context of DFT, Where Does the Nyquist Frequency Sample Belong In a Double Sided Frequency Spectrum (Positive / Negative Side)?

It's convention, they're equivalent:
$$
\exp{\left(j2 \pi \frac{N}{2}n/N \right)}
 = \exp{\left(j2\pi \frac{-N}{2}n/N\right)} \\
\Rightarrow e^{j\pi n} = e^{-j \pi n} \Rightarrow \cos(\pi n) = \cos(-\pi n)=(-1)^n,\ j\sin(\pi n) = j\sin(-\pi n) = 0
$$
MATLAB and Numpy go $[-N/2, ..., N/2-1]$, which is unfortunate for analytic representations (+ freqs only). Note also its value is doubled relative to other bins (but not manually; they correlate this way), so in a sense it's both a negative and positive frequency, so energy's preserved:

You can tell a library's preference by fftshift docs:


Is fft2 in MATLAB unitary? Some differences happen

The standard (conventional) definition of DFT (1D or 2D) is not unitary.

See for example the 1D standard (conventional) DFT pair as:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k } $$

and

$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} k n} $$

DFT is not unitary due to the fact that the forward and backward transforms are not symmetric (the scale $\frac{1}{N}$)

Hence MATLAB's fft() and fft2() functions will not provide unitary transforms; i.e., they do not preserve energy.

However, you can define unitary version of DFT by distributing the scale as:

$$ X[k] = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k } $$

and

$$ x[n] = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X[k] e^{j \frac{2\pi}{N} k n} $$

Now this pair is unitary.


Ways to parallelize Sifting in EMD (Empirical Mode Decomposition)

Is it possible to parallelize the Sifting process?

No.

The "problem" here is not the recursion itself but the fact that the data changes at each recursion step.

It is possible to parallelise divide and conquer type of algorithms that imply hierarchy as is the case of the split radix FFT or Quicksort for example. But in that case, the input sequence of values, let's call it $x[n]$ does not change with every step of recursion. Therefore, this ends up being an "exercise" in managing subsequences and the results of function application. Consequently, it is possible to resolve recursion by storing which indices are supposed to be processed at each step of recursion and storing them on a list, effectively transforming recursion to iteration.

Contrast this with empirical mode decomposition where the vector that is sent to the lower (or higher) levels of recursion is not $x[n]$ but $x[n]-q[n]$ where $q[n]$ is the mean waveform that is formed by combining the maxima and minima interpolated wavefoms.

In this way, the maxima and minima of the second level depend on the subtraction of the previous level. It is not simply a matter of finding the "first order" maxima and minima, then the "second order" maxima and minima, then the third order and so on and then recombining them at different levels of recursion because the first level subtraction has the potential to change which $x[n]$ become local maxima (or minima) at a subsequent step.

Therefore, empirical mode decomposition has an inherent sequential part that cannot be parallelised.

Hope this helps.


Low-pass filtering is based on the assumption that "natural" images have more energy in the low-frequency coefficients than in the high frequency coefficients ; while noise will affect equally all coefficients. Thus, removing the high frequency coefficients will comparatively eliminate more noise than signal. The problem is that there are "legit" high-frequency coefficients in images, such as edges. Denoising by low-pass filtering will attenuate those coefficients and cause edge blurring.

The property used in magnitude thresholding is sparsity - "natural" images are very likely to have only a small set of high non-zero coefficients in the frequency domain. Adding uniform (and independent from pixel to pixel) noise is equivalent to adding a small random value to all the frequency coefficients. The result is that all the coefficients that were 0 in the original image now have a small value ; while the coefficients that were high in the original image are comparatively untouched. Thresholding the small magnitude coefficients will cancel the contribution on noise on those coefficients - though it will not recover the effect of noise on the high coefficients.

If you replace edges by transients/attacks, the same thing applies to sounds by the way.


Many people have asked themselves this question. Unfortunately, there is no exact solution, just approximations. One relevant paper is Practical Design Rules for Optimum Finite Impulse Response Low-Pass Digital Filters by Herrmann, Rabiner and Chan, where they basically fit a formula to their data obtained from lots of designs. A simplified formula is [1]

$$M=\frac{-10\log_{10}(\delta_1\delta_2)-13}{2.324\;\Delta\omega}$$

where $M$ is the filter order, and $\Delta\omega=\omega_s-\omega_p$ is the width of the transition band.

Also the Matlab function firpmord.m uses some formula to estimate the necessary filter order, but they also say that the required order is often underestimated.

[1] Discrete-Time Signal Processing, Oppenheim, Schafer, Buck (2nd ed.), p. 502


why derivative and logarithm operators amplify the noise effects?

The derivative operator is linear and time-invariant and can be described by its frequency response

$$X(j\omega)=j\omega\tag{1}$$

From (1) it is obvious that its gain increases linearly with frequency, and high frequencies (where the SNR is usually bad) are greatly amplified (as already pointed out in endolith's comment).

The logarithm is a memoryless nonlinearity and it compresses the input signal. This operation is used in many signal processing algorithms, such as feature extraction for speech recognition. In these algorithms this operation is essential to the performance and it definitely does not deteriorate the signal. What is important, however, is to deal with the singularity of the logarithm as its argument approaches zero because: $$\lim_{x\rightarrow 0+}\log(x)=-\infty$$
If nothing is done about it, the logarithm will produce very high negative peaks during periods where the signal is almost zero (i.e. when there's only noise). This problem is easily handled by defining a threshold, either on the value of the logarithm or on the input signal. Note that the logarithm can only be applied to non-negative signals (usually magnitudes or squared magnitudes of frequency domain data).


Implementation of a filter in a software application

If your filter is not of insane order, than double precision (64bit on most machines) floating point will totally suffice, and you really shouldn't worry too much. I agree with Arnfinn here.

I don't 100% agree on the optimization recommendation – in fact, if you want my two cents on general purpose CPU-based software implementations of filters,
then simply don't use IIRs; they might be shorter than FIRs with the same flatness/transition band steepness/suppression, but the fact that they are recursive makes them a whole lot harder to execute fast, not to mention that shortcuts like polyphased decomposition don't generally work, and thus, you can't save as much CPU when e.g. building a decimating anti-aliasing filter.

My honest advise is:

Don't worry too much about CPU performance. If you're not trying to do this in JavaScript, and you got but a couple of kilosamples to second to worry about, a simple filter implementation will do.
If you're worried about numerical accuracy, don't worry more than at the level that's sensible: unless your filter is multiple thousand taps long or your taps are extremely bad-behaved, a FIR won't really hurt much. If you're worried about stability, simply don't use an IIR.
Don't reinvent the wheel if you don't have to. There's excellent optimized DSP libraries out there, and as long as the programming language you use supports native libraries with a C calling convention ABI (it probably does), you can use most of them. Just let whatever tool you chose provide the taps for your IIR or FIR, and that library will do its job. I personally like libVOLK a lot, but it concentrates on FIRs, for the reasons given above.

Generally: You're asking about numerical accuracy without indicating which boundaries you need to work in and what filter orders you're talking about. This might indicate that you're not really having any problem at all, because otherwise, you might have noticed already that these are the critical aspects here :)


How to find a 'square root filter' such that $d*\overline{d}=f$ given $f$?

If you can design an FIR filter with arbitrary magnitude response $F(\omega)$, then you can easily design the filters that you're looking for. Just break the design into two pieces: magnitude and phase.

If the desired magnitude response of the cascade is $|F(\omega)| = |D(\omega)||D^*(\omega)| = |D(\omega)|^2$, then, as was pointed out already, you can simply take the square root of the desired overall magnitude response to yield the magnitude response of your square-root filter $D(\omega$).

With the magnitude response defined, then it's up to you to pick what kind of phase response you want to give $D(\omega)$. The phase responses of $D(\omega)$ and $D^*(\omega)$ will be conjugates of one another, so they will cancel when the two square-root filters are cascaded. You would typically pick the phase response to aid in implementation or to give you some other nice property along the way.

One common choice, as you pointed out, is using linear-phase FIR filters. So, you can pick a filter length, set the filter group delay to $\frac{N}{2}$ (where $N$ is the filter order), and use the resulting phase delay curve as the desired phase response in your square-root filter design. You should end up with a filter with a symmetric impulse response (which can provide implementation advantages due to the redundancy in the structure of the taps).

One other thing I should note is that the question as posed will yield a non-causal structure in one of the subfilters. Take the case where $D(\omega)$ has linear phase (so a constant group delay). Its conjugate will have the opposite phase response, so if $D(\omega)$ corresponds to a delay of $N$ samples, then $D^*(\omega)$ will have a phase response that corresponds to an advance of $N$ samples, so it will be non-causal.

What the OP was probably looking for was a structure where $D^2(\omega) = F(\omega)$ instead. This would be similar to the root-raised cosine filter structure used in digital communications, for example.


Since this is homework/a learning assignment, I think a hint is worth more than the (admittedly done on this site a couple of times) complete proof:

You forgot to also read the context: That's not true for just any noise and any signal, it's only true for noise that's white and a signal that fits perfectly into the filter shape.

Think about what "white" has to do with a moving average. Hope this helps.


Is it better to use sampling rates which are a power of 2? 32-64-128-256-512Hz? instead of using for example 30- 60-100-250-500Hz?

There's a lot of things that influence your choice of sampling rate.

First of all, of course, the Nyquist theorem, which says you need to sample at more than twice of the signal's single-sided bandwidth $B$ (for real signals)¹.

Then, you often choose a sampling rate that is useful for your system. For example, a lot of 3G/4G/5G systems use sampling rates of 30.72 MHz (or multiples thereof). That's not an elegant multiple of anything (as far as I know) – it just happens to be what the symbol rates in the standard suggests, so it's very useful if you're already sampling at rate that is rational relative to what you'll be using in the end.

That's a generally important rule: Pick a sampling rate that makes sense for what you'll be doing to the signal afterwards.

Then, especially in audio and lower-frequency sampling, you choose sampling rates that are multiples of other rates, so that it's easy to keep multiple subsystems in synchronization. For example, your CD audio sampling rate is 44100 Hz – first of all, that's an integer multiple of 100 Hz, and that's good, because it means it's easy to keep not only the second counter on your CD player correct by simply counting the samples and dividing by 44100, but also things like studio hardware, which tends to have higher "human-readable" time resolution. But, why the strange rate 44100? Well, that rate happens to be $44100=2\cdot2\cdot3\cdot3\cdot5\cdot5\cdot7\cdot7$. Back in the 1980s, when that rate came up, it was harder than it is now to derive good clocks from each other, unless they had an integer factor between them. So, a clock with "all" the small integer factors was chosen. A curious choice!

Another reason to pick a specific sampling rate is noise alignment, or intentional non-alignment. While it's not a sampling rate itself, the NTSC frame rate of 29.97 frames a second is still something that you find in video systems today, so you find a lot of video capturing ADCs that run at a high multiple of 29.97 Hz. And that frequency was chosen so that the zeros of the color information shaping filter didn't interfere with the brightness info.

In the very same motivation, you'll find noise-sensitive systems where the ADC rate is as far as possible from being rationally related to the other dominant digital clocks in the system – so that the noise from these clocks distributes more evenly to the digitized signal.

In other systems, you painstakingly try to align all these clocks – so that clock noise ends up being narrow interferers in the digital signal, and thus can be digitally filtered out using notch filters.

Of course, you also need to be able to get any rate you'd want to use – you can't just by a 876.543 kHz oscillator, even if you wanted to, these don't exist commercially. Sure, you can synthesize one with either a lot of jitter or a lot of effort, but in the end, complexity and cost restraints are relevant for almost any system design.

As you can see, a lot of conflicting reasons!

¹ that's not necessarily twice the maximum frequency in the signal! A lot of systems, especially superheterodyne systems, i.e. systems that have a strictly bandslimited signal on an intermediate frequency (IF) use subsampling, i.e. the sample rate is $> 2\cdot B$, but $\ll 2\cdot f_\text{max}$. The aliasing that occurs doesn't distort the information content, but aliases the signal down from the IF to the Nyquist band. It's a very useful, very common technique.


What is the ROC for this discrete signal:

HINT:

Remember that the ROC is the region in the $z$-plane for which the series

$$X(z)=\sum_{n=-\infty}^{\infty}x[n]z^{-n}\tag{1}$$

converges. Since you've found $X(z)$ you also know that $x[n]=4\delta[n-2]$, i.e. there's only one value of $n$ for which $x[n]$ is not equal to zero. What does that mean for the convergence of $(1)$?


Ideas on Matrix Factorization / Transformations for $ {L}_{1} $ Minimization

Since $ \epsilon $ is a parameter you need to set why not trade it with another parameter you need to set to create an easily solvable problem (Relaxation of the Problem)?

You can transform the problem into the following form ($ {L}_{1} $ Regularized Least Squares):

$$ \arg \min_{x} \frac{1}{2} \left\| A x - z \right\|^{2} + \lambda \left\| x \right\|_{1} $$

There are many solvers for that even for large sparse matrices.
For any $ \epsilon $ from the original problem there is $ \lambda $ such that the solutions are the same.
Since $ \epsilon $ it to be arbitrarily set in the original problem you can trade it for $ \lambda $.

This can be solved with many solvers easily without the need for Matrix Inversion.


Continuum Fitting with Python using specutils

At the end I settled on using a median filter from scipy (medfilt) and with a high enough kernel it removed the spikes completely.
The problem on my specutils implementation was that I needed to exclude the peaks in the spectrum window by doing this, althought I found this too hard to implement on 100-120 data measurement I had:
...
from specutils import SpectralRegion

spectrum = Spectrum1D(flux=y*u.Jy, spectral_axis=x*u.um)
g1_fit = fit_generic_continuum(spectrum, exclude_regions=[SpectralRegion(2 * u.um, 6 * u.um), SpectralRegion(58 * u.um, 65 * u.um)])
y_fit = g1_fit(x*u.um)


how I can use the FFT data to compute the convolution?

You can pad zeros in frequency domain. But you need to take care about Nyquist bin when N is even.

A3 = [A2(1:N/2) A2(N/2+1)/2 zeros(1, 2*N-length(A2)-1) A2(N/2+1)/2 A2(N/2+2:end)];

Here Nyquist bin is A2(N/2+1)/2 and zeros are padded in the middle of the spectrum

For N-odd

A3 = [A2(1:(N-1)/2+1) zeros(1, 2*N-length(A2)) A2((N-1)/2+2:end)];

This is only for the case if you need to interpolate in time-domain. Almost irrelevant to the problem you have.

Posted the solution to you problem below.


I'm assuming since you didn't say that you're working with audio. But my answer generally applies if you're working with images/video as well. Unless the background you're trying to subtract is exactly the same, I doubt subtracting a phase-inverted version of the signal will work. It sounds like what you're trying to perform is spectral subtraction. This is a common technique where you simply subtract the average FFT magnitude spectrum from your signal. See this paper on spectral subtraction. If working with images/video, I assume that you could do this in 2 dimensions as well and the math would be similar.


What is the relationship between frequency vector and the peak frequency in FFT?

Although I think Jason's answer contains all the information you requested, maybe what you're looking for is a little less theoretical.

Let's assume you have a signal of $$x(t)=cos(2\pi t f_x)$$
where $t$ is the time in seconds and its frequency $f_x=20$Hz. Let's assume further you sample this signal at $f_s=1000$Hz over $T=3$ seconds, i.e. you have $K=3000$ samples $x[k]$, where $k$ is the integer sample index of the signal $x(t)$.

The sampling theorem basically says that your sampling frequency $f_s$ must be at least twice the maximum signal frequency (which is the Nyquist frequency $f_s/2$) in order to reproduce it without aliasing. For our example: $f_s=1000>2f_x=40$.

You choose the "number of points" $N$ for the FFT (you called it NFFT) yourself, depending on how many different frequency samples you need to distinguish. The maximum frequency resolution you get by using all available samples you have for calculating the FFT, which - depending on the $f_s$ and $T$ of your signal - might not be feasible.

For our example, using all 3000 samples $x[k]$ is not a problem on a PC running Matlab, so let's assume we do that get $N=3000$ complex frequency samples $X[n]$ where $n$ is the integer frequency index.

Now to your question where in $X[n]$ to find a specific frequency $f$ you're interested in if $N=K$: $$n = f T$$

Therefore, the Nyquist frequency $f_s/2$ bin is at $n_s = T f_s/2 = 1500$, right in the middle of the FFT result vector. If you plot the amplitude spectrum $abs(X)/N$, you'll notice that there are peaks at $n_x = T f_x = 60$ and $n_x' = 2940 = N - 60$, each of magnitude $1/2$. This is because the spectrum is "mirrored" around $f_s$ and the signal amplitude of 1 is divided between both peaks $n_x$ and $n_x'$.

For Matlab indices, I think you'll have to add 1 to the $n$ mentioned above because Matlab's vector indices start at 1 instead of 0.

Note that if you have frequencies $f$ which are not a multiple of your sampling frequency $f_s$, the amplitude will spread to neighbouring bins and you'll have to round $n$ to the nearest integer.

See also Matlab's FFT documentation, which has a nice example.


What Is the Difference Between PCA and Karhunen Loeve (KL Transform)?

For discrete data both are the same - Finding set of orthogonal directions which maximizes the Variance (Energy) of data along them. Sometimes those are called the natural axis of the.
Since we're dealing with variance it is only natural both are calculated from the covariance matrix of data.

You may encounter places where KL might be even defined on the correlation matrix.

Anyhow, some nuances when dealing with them:

Usually when dealing with Continuous Function people use the term KL. PCA is preserved to discrete data (Though Mathematicians might work with the PCA in all cases as they have SVD for continuous data as well).
PCA is used for any discrete data where KL was specifically formed with random variables in mind.

In practice, when one has discrete data and tries to find the directions which maximizes the energy of the data (Dimensionality reduction, compression, etc...) both are the same.


Resampling: how many samples to zero-stuff or down sample?

When you up sample by 160, if you wanted to filter twice, you would actually be applying an "anti-image" filter.  The next filter for decimation would be the "anti-alias" filter.  But in you case, you could just use the anti alias filter once and accomplish the same.  I think perhaps you need to consider what zero stuffing means, to put so many zeros between samples and decimate is to throw away so many samples inbetween.

Btw, you do compute samples for the zeros when you filter.  If you want to understand better, try plotting the spectrum of each stage, and what would happen if only using the 13.

If I understand your question properly, you would effective just upsampling by 13, filter, then downsample.  You'd be back where you started


texture matching between patches of an image

The mean and standard deviation are two measurements of a distribution. Others you could also use are higher order moments like 'skewness' (how skewed the distribution is) and 'kurtosis' (how 'peaky' the distribution is).

However, what I would try is a histogram of the values for each of the channels. For example, if you used a histogram with 16 bins, you would get a feature vector with 48 values. The histogram bin values can be thought of as a descriptor of the distribution.

To account for changes in brightness, you need to introduce some illumination invariance. Already the LAB colour-space splits the image into luminance (L) and color (AB) channels. To achieve illumination invariance I would normalise the L channel histogram values using the norm of whatever measure you are using (in this case, euclidean distance, it would be the l2 norm). Normalising the A or B channels is probably not necessary. The normalisation accounts for any linear changes in luminance.

Depending on the images, you could also add vectors like:

histogram of oriented gradients (HOG)
histogram of connected components (make B/W then count how many pixels are white in the 8 connected neighbours)
k-means clustering of LAB values to get k centre LAB values which you can compare to other images.

and so on.


python: how to compute the gray level histogram features as mentioned in the paper, and

(1) how to we compute the gray level histogram Hr , Hg , and Hb for each of the red, green, and blue channels.

The gray level histogram for each of the channels is simply the gray level histogram of the Red channel, the Green channel and the Blue channel separately. Instead of converting the image from RGB to grayscale, try to extract its RGB components.

For example, in GNU Octave, a colour image is represented as a $M \times N \times D$ matrix, where $M$ is the height of the image, $N$ is the width of the image and $D$ is the "number of channels".

In a colour RGB image, the number of channels is 3.

This can be viewed as 3 images packed together in one "3D" matrix. To work on the Red channel, all that you have to do is use slicing to obtain just one of those matrices.

In Octave, this would be done via:

image = imread("someimage.png")
red_channel = image(:,:,1)

Similar thinking applies to other languages / computing platforms.

  2) how to compute the the contrast quality, qct mentioned in the paper.

According to the paper, the $q_{ct}$ is the width of the 98% "part" of the combined histogram.

So:

Load the image
Get the histogram of the Red Channel (call it $H_r$)
Get the histogram of the Green Channel (call it $H_g$)
Get the histogram of the Blue Channel (call it $H_b$)
Produce the combined histogram as $H = H_r + H_b + H_g$
Get the 98% of its mass:

Sort $H$ in reverse order of values (largest count comes first)
Get the cumulative sum of the histogram values until you have covered the 98% of the total sum.
Find the minimum and maximum values of the indices of the values that contribute to the sum
Calculate the $q_{ct}$ as the difference between the maximum and minimum of the indices.

Hope this helps.


Let me slightly change notation and denote the two modulation frequencies by $\omega_1$ and $\omega_2$.
You need $|\omega_1-\omega_2|>2\omega_{max}$ to be satisfied, otherwise the shifted spectral of $x(t)$ and $y(t)$ will overlap. If this condition is satisfied, you can retrieve $x(t)$ by demodulation (multiplication with $\cos(\omega_1t)$) and (scaled) low pass filtering, as you suggested:
$$\begin{align}2x(t)\cos^2(\omega_1t)=x(t)\big[1+\cos(2\omega_1t)\big]\Big|_{low pass}=x(t)\end{align}\tag{1}$$


• The N-Ifft points and number of subcarriers are not same all the time. Why we introduce unmodulated or zero padded subcarrier between the modulated or useful subcarrier.

For different reasons. Mainly:

Most OFDM systems are based on direct-conversion receivers. In those receivers, you mix a signal down to complex baseband by multiplying it with a complex tone of the signal's center frequency – but that leads to the signal at exactly 0Hz in the baseband to be unrecoverable, technologically. For that reasons, the DC carrier (or carriers around DC) (DC=direct current, a very Electrical Engineering way of saying "0Hz") is not used
OFDM systems often have to deal with problems of strong Out-Of-Band radiation. By leaving the "edge" carriers unused, you can mitigate that.

  • Should the training sequence be added before ifft operation or after ifft+cp inclusion operation? If both can be possible then what is the best approach between these two.

I assume you want to train your equalizer? There's a lot of approaches to channel equalization with OFDM, and every one has its advantages and disadvantages.

The same applies to adding a timing and frequency synchronization sequence: there's a lot of approaches to doing synchronization for OFDM, and hence, as many trade-offs for advantages and disadvantages to having something in time- and frequency domain. Note that fine frequency estimation is often done with known time-domain signal techniques (because we can't do an FFT overly well if we don't even know where), and coarse frequency correction (ie. correcting an offset of integer multiples of the subcarrier spacing) can then be done afterwards (in time- or frequency domain, whichever is more appropriate for the system.

I personally would recommend Schmidl&Cox' classic paper Robust frequency and timing synchronization for OFDM – it's a great read, and I personally think it's a good paper, with a short, concise intro, a well-explained method, well-condensed mathematical proofs, and useful explanation and derivation of performance parameters. Also, if you're doing OFDM in your thesis, and want to actually have read at least one of the papers you cite, this should probably be the one.

  • For the DAC conversion, raised cosine pulse shaping filters are considered to be the best among other shaping filters. What should be its roll off factors, No. of symbols and No of samples in each symbol (for OFDM). How much should be an oversampling factor for upsampling process.

Generally, as you've learned, Roll-off etc are parameters of your design that you choose to make your design work according to your requirements. So, no matter for which method you'd be asking, "what are the values of the Roll off factor for ?", the answer will always be, "what you design them to be!".

However: OFDM was partially "invented" to eliminate the need for a set of pulse shaping filters. So your question regarding a raised cosine doesn't make much sense. If you pulse shape the time domain OFDM signal, you destroy the "O" in OFDM! You mustn't do that.

  • Is Low Pass Filter should be used at the transmitter side too what is the purpose of using it at the transmitter side ?

I think you're sitting in front of some specific example, and think it is "all possible OFDM transceivers". I can't tell you what your example is – I don't know it. Maybe that low-pass filter is just part of the direct conversion mixer that converts the OFDM baseband signal to a passband signal?!


As mentioned in a comment, what you describe looks like a timing recovery algorithm for digital communications.
I can recommend two sources, an article and a chapter in a book, which give a nice introduction and overview of the topic in general:

L. R. Litwin, «Matched filtering and timing recovery in digital receivers», RF Design, pp. 32–48, Sep. 2001.

M. Rice, «Chapter 8. Symbol Timing Synchronization», in Digital Communications: a Discrete-Time Approach, Upper Saddle River, N.J: Pearson/Prentice Hall, 2009.


Passband attenuation refers to the maximum signal loss within the passband of the filter which should ideally be zero; i.e., pass band shall not attenuate the signal it passes; it's a unit gain region.
The term ripple is about of alternating gain through the passband for certain filter types yet it's still used even if the passband is not alternating such as a Butterworth filter. In that case the frequency at which the gain reduces to -3dB is the cutoff frequency for a lowpass filter. Or similarly high and low cutoff frequencies for a bandpass filter.


Why is an RC circuit a first order system?

The order of system gives the number of memory storing elements in the circuit.

What you have written above is a KVL equation of the circuit. But when you determine the order of the system you need to find the transfer function of the system. Laplace domain is normally used. The highest power of $s$ in  the denominator in a strictly proper transfer function gives the order of the system.

So in the $s$-domain,say you take the output across the capacitor,

$\dfrac{V_o(s)}{V_i(s)} = \dfrac{\frac{1}{sC}}{R+\frac{1}{sC}}$

simplifying we get

$\dfrac{V_o(s)}{V_i(s)}= H(s)=\dfrac{\frac{1}{RC}}{s+\frac{1}{RC}}$

So the highest power of $s$ is given by 1. So the order of the system is 1.


The essential bandwidth of a rectangular pulse

First note that $$\mathcal{F}\{\mathrm{rect}(\frac{t}{T})\}=T\mathrm{sinc}(Tf)=\frac{\sin(\pi Tf)}{\pi f}$$
Then you have correctly calculated the pulse energy as $$E_g=\int_{-\infty}^{+\infty}g^2(t)dt=T$$
Now we want to find $B$ such that
$$\int_{-B}^{B}T^2\mathrm{sinc}^2(Tf)df=0.9T$$
or since the function is even,
$$\int_{0}^{B}T\mathrm{sinc}^2(Tf)df=0.45$$ which by assuming $$Tf=x\Rightarrow df=dx/T$$ is
$$\int_{0}^{TB}\mathrm{sinc}^2(x)dx=0.45$$
Using the following code we get the essential bandwidth $$\bbox[0.5em,#efe,border:0.1em groove navy]{\
B\approx\frac{0.8455}{T}\
}$$

J = 0.45; %target value
acceptable_error = 1e-4;
dx = 1e-5;
xmax = 1; %initial guess
max_iterations = 100;
iteration = 0;
while iteration < max_iterations;
    x = 0:dx:xmax;
    y=sinc(x).^2;
    I = trapz(x,y);
    if abs(I-J) <= acceptable_error
        break;
    end
    if I < J
        xmax = 3/2*xmax;
    else
        xmax = xmax/2;
    end
    iteration = iteration+1;
end
fprintf('absolute error = %f\n',abs(I-J));
fprintf('solution = %f\n',xmax)


How to distinguish between recorded and a live vocal signals?

On a smartphone, flat out no. The microphones are so poor that you're lucky to get any approximation of the actual audio.

On high-quality hardware, perhaps. Barely. If you're lucky. Check for MP3 compression artifacts. If the source is not compressed with a lossy codec, you're out of luck.


OpenCV Motion Analysis: Optical Flow vs Motion History

Farneback is a dense optical flow algorithm. This means, it outputs the flow vectors per each pixel (which can be tracked). LK is a sparse variant, tracking only certain feature points.

The motion segmentation algorithms generally target estimation of foreground or background and doesn't allow establishing correspondences. Please refer to this post. You will find more details in my response.


Converting from RGB to LAB Color Space to Reduce Colors

You should read about LAB Color Space.

Its main advantage on RGB is by having 1 channel dedicated to the Luminosity of the Image and 2 others dedicated to the color information.

This means if you need to process only Luminosity / Color you can do that using LAB without having any effect on Color / Luminosity.


(IANAL...)

If you only want detectors:

Harris is probably OK.
According to http://users.fmrib.ox.ac.uk/~steve/susan/, SUSAN is out of patent now.
I've not seen any claims that FAST is patented.

Descriptors are harder...

Histograms of Oriented Gradients might be worth considering - again I've not seen any claims of patent on the original form.


Is the frequency spectrum dependant on frequency, or on imaginary angular frequency?

Basicaly it's just a convention and the Fourier transform is a function of $\omega$ .

The imaginary argument $j\omega$ can be understood in different perspectives.
In one convention, for the class of signals $x(t)$ which are absolutely summable (i.e., stable) , the continuous time Fourier Transform (CTFT) $X(j\omega)$ is considered to be derived from the existing Laplace transform $X(s)$, with $s$ being a complex number, by replacing $s$ with $j\omega$, which amounts to evaluating the Laplace transform on the imaginary axis ( where $s =j\omega$ )  of the complex plane. Hence in this interpretation, the argument $j\omega$ retains the conceptual and mathematical link between the Laplace transform $X(s)$ and the CTFT $X(j\omega)$ where we say $X(j\omega) = X(s)|_{s=j\omega}$

A similar argumentation also holds for the discrete time Fourier transform (DTFT) $X(e^{j\omega})$ as such: For the class of signals $x[n]$ which are absolutely summable (i.e., stable), the discrete time Fourier Transform $X(e^{j\omega})$ is considered to be derived from the existing $Z$-transform $X(z)$, where $z$ is a complex number, by replacing $z$ with $e^{j\omega}$, which amounts to evaluating the $Z$-transform on the unit circle ( where $z = e^{j\omega}$)  of the complex plane. Hence in this interpretation, the argument $e^{j\omega}$ retains the conceptual and mathematical link between the $Z$-Transform $X(z)$ and the DTFT $X(e^{j\omega})$ where we say $X(e^{j\omega}) = X(z)|_{z=e^{j\omega}}$

However, note that, there exist certain very important class of signals (such as the pure sine wave or the impulse response of an ideal Low-Pass filter) which are not absolutely summable (unstable) therefore whose formal Fourier transforms (as obtained by replacing $s$ or $z$ by the corresponding $j\omega$ or $e^{j\omega}$) do not exist. Those Fourier transforms are obtained from another perspective by utilizing the impulse function $\delta(\omega)$ or discontinuities in the frequency domain. In such a case the argument is misleading and there exist no coresponding $X(s)$ or $X(z)$ which you can obtain by replacing $j\omega$ or $e^{j\omega}$ with $s$ or $z$ in the Fourier transform $X(j\omega)$ or $X(e^{j\omega})$, hence in that case it's just a convention.


Unexpected behaviour when using lfilter with initial filter delay multiple times

The problem is not in the filtering process but already at the design stage. Your specifications are very difficult to realize because your desired band is at very low frequencies. With these specs, the butter routine returns an unstable filter (i.e., with two pole pairs outside the unit circle of the complex plane).

One thing you could try is reduce the filter order. Always check the maximum pole radius by computing the (magnitude of the) roots of the denominator polynomial. It could be that the design turns out to be OK for a lower filter order. Of course, your filter will be less steep.

Another, better approach would be to downsample your signal before filtering, so your pass band is not anymore in such a low frequency range (compared to the sampling frequency). This will avoid numerical problems in the filter design routine.


How do I prove that delta - sinc function is the same as an (-1)^n times the sinc

Hint : Use the fact that, $\sin(-\frac{n\pi}{2}) = (-1)^n \sin(\frac{n\pi}{2}), \forall n \in \mathbb Z$


If the frequency response of a system is time-independent, is this system LTI?

LTI means Linear Time-Invariant systems. The system has to satisfy two conditions.
(1) Linear and (2) Time-Invariant.
(1) Linear means, if the response of the system due to load Px and Py is Rx and Ry respectively, then for the load (Px+Py), the response of the system will be (Rx+Ry).
(2) Time-Invariant means, the parameters of the system does not vary with time.

If the frequency response of the system is Time-Independent, then it suggests that the parameters of the systems are not changing over time, i.e., the system is time-invariant.

But it does not say anything about nonlinearity. Hence, in a nutshell, if the frequency response of a system is time-independent, this system is NOT an LTI?


Convolution of discrete signals using convolution sum

Well, you know that $u[n]=0$ for all $n<0$. So from the two terms in your sum, all terms will be zero where either $k+2 < 0$ (i.e., $k < -2$) or $n-2-k<0$ (i.e., $k>n-2$). From this your sum will need to run from $-2$ to $n-2$. Can you take it from here?


Deriving the CRLB for an active system with non-square shaped bandwidth

Found the solution, I believe.

It's a dumb mistake. I'm calculating the SNR as $S_0 / N_0$ in one case and as $S/N$ in the other.

Since the noise is square shaped, and the signal is triangle shaped, the first is obviously wrong.
Correcting to $\frac{S_0 B/2}{N_0 B}$ gives the factor that was missing from one to the other.


I thought about it a bit more (power of posting), and decided the best way to scale $x$ and $\dot{x}$ is by the standard deviation, which is being tracked within the covariance matrix.

errX = (x_n - x_n-1)/sqrt(var(x));
errXdot = (xdot_n - xdot_n-1)/sqrt(var(xdot));
err2D = sqrt(errX^2 + errXdot^2);


There are many options to test a wavelet transform, or compare them.
You can first compare them on their most basic properties. One can check whether it is invertible, orthogonal, biorthogonal, critical or redundant (and the oversampling ratio), separable or not, (close-to) invariant to transformations of the signal (symmetry, rotation, scaling, shear). They take weak care about data. Because those properties are often tested against "all data", or random enough (like several noise realizations).
And at the end of the processing workflow (denoising, compression, restoration), you can compare the outcome of some objective or subjective measures (SNR, PSNR, SSIM, other norm-related metrics).
In between, one can estimate vector features on wavelet coefficients. Such measures can be computed and evaluated  via standard statistics (mean, standard deviation, median, median absolute deviation) and across data models. But you can also think of higher order moments: absolute, centered, ratios of moments. From now on, I will call coefficients from the same scale "subband" (aka voices for continuous wavelet transforms), and use the discrete wavelet transform formalism. Beware: subbands can be subsampled in different manners, and scaled accordingly.
Typically, one computes a vector made of statistics of each subband, or a gathering of subbands. As wavelet subbands are normally zero-mean, coefficients are typically norms (non-centered moments) of subband coefficients. Suppose a three-level decomposition, with $A_3$, $D_3$,  $D_2$ and $D_1$ approximation and detail coefficients. A feature vector could be made of the norms of $D_3$,  $D_2$ and $D_1$, and may something for $A_3$ (but for standard wavelets, the approximations could be close.
Sparsity is another interesting feature: how well are wavelet coefficients concentrated with respect to the original signal? A 2009 paper, Comparing Measures of Sparsity, provided axioms to compare statistical distributions, among them: Gini indices, norm ratios, kurtosis, etc.
As they have been recently considered for quality metrics, stopping criteria or robust sparsity estimators, I would like to pinpoint two of our recent adventures into using regularized  ratios of norms ($\ell_1/\ell_2$ in Euclid in a Taxicab: Sparse Blind Deconvolution with Smoothed l1/l2 Regularization) or SPOQ ℓp-Over-ℓq Regularization for Sparse Signal Recovery applied to Mass Spectrometry with its Matlab toolbox.
Norm or quasi-norm/norm ratios can be quite efficient while dealing with sparsifying transformations.


Calculate sum of absolute difference per pixel

I think there are multiple things to consider here:

1) Where are you dividing the image into 8x8 blocks? If you don't divide like this, then normalization by 64 is meaningless, and you are just computing SAD over all pixels. Maybe I'm seeing something wrong? You should do the normalized SAD in those small windows and then sum and normalize. And since 3600 is not a normalized value, I don't think that he normalizes.

2) Have you considered that your function for loading the video or the codec you use is actually giving you slightly different values than the one he uses? If your input is different, you will obviously get different results.


How to initialize an IIR filter in Matlab

Not easy. The Matlab filter() function implements a general purpose IIR filter and allows you to pass in an initial state. However, they don't publish the exact definition of the state and the algorithm deployed. It seems to be close to a transposed form II filter but with some subtle numerical differences, so the state doesn't quite match.

In any case it would be better to break this down into second order sections using sosfilter(). However that doesn't allow you to set an initial state. You may have to write your own filter function for that.

YOu still can use filter() to some extent. The default is zero initialization. You can also do step initialization by simply running a constant input for a while and collecting the state as an output argument from the filter, like this:

[b a] = butter(4,1000*2/44100);
% create state for a constant input of "1"
[~,state1] = filter(b,a,ones(1000,1));
% now filter for a constant of 2 but with a state for "1" input
ystep = filter(b,a,2*ones(100,1),state1);
plot(ystep);

Not sure about projection initialization, though.


and how to avoid it?

Address the NaN in your input data. Don't "fix", "paper over it" or "replace by 0". Find the root cause for the NaN, understand what's happening and take meaningful corrective action.

NaN means your input data is bad. Doing anything with bad data is pointless since your output will be bad. Doing cosmetic adjustment just so the code doesn't blow up is equally pointless. You need to understand why the data is bad and address the root cause


What is the state of the art for vocal elimination/reduction?

Some recent ideas (2-3 years old, I haven't followed recently):

Design an elaborate generative model of the spectrogram of music signals, combining a non-negative mixture of a small set of background frames (models the background music) + a harmonic comb filtered by a non-negative mixture of smooth basis functions (the source filter model of the lead instrument), with constraints on the pitch transitions. Estimate model parameters from the audio clip. Use Wiener filtering to resynthesize the background music/vocals from the estimated spectra. See A musically motivated mid-level representation for pitch estimation and musical audio source separation by Durrieu et al. Such Bayesian approaches can be made more complex ad nauseam, by integrating prior knowledge about note transitions/durations - at the cost of increasingly hairy inference procedures. Demo time!
Techniques inspired by background subtraction in video analysis. The idea is that the accompaniment (especially in "bland" pop music) is invariant by a translation of a few bars (the duration of the "loop" or "riff") - while the vocals are not. Once the structure of the piece of music is known and can be described in terms of repeated groups of bars (a well-studied topic in itself), one can take the median of translated versions of the spectrogram to suppress the vocals. See REpeating Pattern Extraction Technique (REPET): A Simple Method for Music/Voice Separation by Rafil & Pardo. Demo time!

In the commercial space, you can download a trial for ADX Trax.


Fair performance comparison betweem LMS & NLMS

For fair comparison of one algorithm to another, the value of step size does not need to be same. You can adjust the step sizes of both algorithm so that the mean-squared-error learning curves base floor gets same and in this way you will be able to differentiate the performance of algorithm. For base floor I mean the value at which the mean squared error curves coincides with each other at the end of iterations or data samples. You can adjust the step size by hit and trial method.


Number of samples in one period of subcarrier after IFFT

If you'd like to think analog, an OFDM signal can be written as a sum of weighted complex sinusoidals:
$$
x(t)=\sum_{k=0}^{N-1}X_k \exp{\left( \mathrm j 2\pi\frac{kf_\mathrm{s}}{N}t \tag{1}\right)},
$$
where $N$ is the number of subcarriers, $f_\mathrm{s}$ is the sampling frequency and $X_k$ denotes the subcarrier values. For a digital implementation, (1) is sampled at $t=n/f_\mathrm{s}$:
$$
x(n/f_\mathrm{s})=\sum_{k=0}^{N-1}X_k \exp{\left( \mathrm j 2\pi\frac{kn}{N} \tag{2}\right)},
$$
which corresponds to the inverse discrete Fourier transform (IDFT), implemented by an IFFT. The IFFT has $N$ input values and $N$ output values. That's ok, because the IDFT is N-periodic and thus any additional sample would be redundant. Consequently, every subcarrier in (1) is represented by $N$ samples, regardless of its number of periods during the OFDM symbol length.

If you try to imagine the digital OFDM time domain signal as a superposition of $N$ sampled sinusoidals, the figure you posted has two problems, in my opinion:

The first subcarrier in the image obviously has frequency $f_\mathrm{s}/N$ and the fourth subcarrier has frequency $f_\mathrm{s}$. That's wrong. From (1) you can see that the first subcarrier should have frequency 0 and the fourth should have frequency $(N-1)f_\mathrm{s}/N$. You can reorder the subcarriers (e.g. from negative to positive) but the DC subcarrier must be included and there must not exist a subcarrier with frequency $f_\mathrm{s}$. (One could argue that subcarriers with frequency 0 and $f_\mathrm{s}$ are actually the same but I think it's a poor choice for an illustration)
The image only shows the real (or imaginary) part of the sinusoidals. That could give you the impression that some aliasing is going on when sampling. I will explain in a second what I mean.

Please consider the following figure that I have created instead

It contains the subcarriers from $k=0$ to $k=3$ (lines) and their samples (dots). It also shows the real part (blue) and imaginary part (red). Now consider the red line in the first and the third subfigure: their samples are exactly the same. However, the blue (imaginary) part is different. That is, $N=4$ samples are sufficient here to represent all subcarriers.

Here is the Matlab code I used to create the above figure:

N = 4; % num subcarriers
os = 128; % oversampling
t = linspace(0, 1-1/(N*os), N*os); % time
n = 0:N-1; % discrete time

figure;
cnt = 1; % helper
for k = 0:N-1;
    subcarrier = exp(1i * 2*pi*k*t);

    subplot(N,1,cnt);
    plot( t, real(subcarrier), 'b-');
    hold on;
    plot( t(n*os+1), real(subcarrier(n*os+1)), 'b.', 'MarkerSize', 20);

    plot( t, imag(subcarrier), 'r--');
    hold on;
    plot( t(n*os+1), imag(subcarrier(n*os+1)), 'r.', 'MarkerSize', 20);

    if(k ~= N-1)
        set(gca, 'XTick', []);
    else
        xlabel('t/T');
    end
    ylim([-1 1]);
    title( sprintf('k=%d', k));

    cnt = cnt + 1;
end


Unsupervised learning algorithms to detect anomaly in waves

Try Nonparametric Statistics.
You may divide your random trajectories into parts and compare their statistic characteristics (mean and variance). Smth like this (but not exactly):

http://www.cscjournals.org/manuscript/Journals/SPIJ/volume7/Issue1/SPIJ-210.pdf


Multiplication with $s$ in the Laplace transform domain equals differentiation in the time domain. In the discrete-time domain we can approximate differentiation by the equation

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

where $T$ is the sampling interval.

In the Z-transform domain, Eq. $(1)$ becomes

$$Y(z)=X(z)\frac{z-1}{T}\tag{2}$$

I.e., the transfer function

$$H(z)=\frac{z-1}{T}\tag{3}$$

approximates differentiation, and replacing $s$ in a continuous-time transfer function by $H(z)$ is thus a way (usually not the best one) to approximate a continuous-time system by a discrete-time system.

One of the problems with the forward Euler method is that transforming a stable continuous-time system could result in an unstable discrete-time system.