How to Generate a vector map for histogram equalization

Your code is halfway correct; it creates the cumulative distribution function (CDF), but the map isnt quite there yet.

First you get the histogram using imhist(), but we want to convert this to probabilities (the likelihood that each greyscale value occurs at a pixel). mapVec gives us integer counts, so just divide by the number of pixels in the image.

mapVec=imhist(Im);
mapVec=mapVec / prod(size(Im));
EqVecMap = mapVec;

The next part calculating the CDF is correct (except for mapVec(i), I changed that to mapVec(j))

for i = 1:size(mapVec)
    EqVecMap(i)=0;
    for j = 1:i
        EqVecMap(i)=EqVecMap(i)+mapVec(j);
    end
end

But we're not done yet. EqVecMap is a CDF, so it's a double vector with max value is 1. However, our pixel values are INTEGERS with max value 255. So just

mapVec=round(EqVecMap*255);

mapVec is now a mapping function from old greyscale value (index 'i') to a new value (mapVec(i)).

** I just saw havakok's answer (but I cant comment bc reputation), and I think it does basically the same thing. Only difference is you need to multiply by 255, not 256, since our max greylevel is 255 (the range includes 0, remember).


Accessing Maximum Value from a Singular Value Decomposed Matrix

The SVD Decompose the image into the (One way to look at it) many matrices.
For instance, given an Image $ I $ its SVD is given by:

$$ I = U S {V}^{T} = \sum_{i=1}^{\textrm{rank}(I)} {s}_{i} {u}_{i} {v}_{i}^{T} $$

As you can see, the image is decomposed into linear combination of matrices (Given by $ {u}_{i} {v}_{i}^{T} $) where the weight of the $ i $ -th matrix is given by $ {s}_{i} $ which is the $ i $-th Singular Value (Usually they are organized in descending order).

The dominant features are usually goes with the large Singular Values (Moreover, the smallest one usually describe noise).
I guess methods want to focus on the large dominant features and avoid dealing with noise hence it makes sense to deal only with the dominant matrix in this decomposition.

By the way, I'd experiment with trying to work on the next ones as well.
In addition, on the recompositon process I'd neglect the last few to have "Simple Denoising".

One should be aware that the SVD of the whole image might be time and memory consuming when dealing with large images.

Moreover, it makes sense to stack few image patches into large matrix and then apply the SVD.
This might extract what is common to all patches (Linear Dimensionality Reduction).


IIR filtering vs IFFT filtering by zero binning and making bandpass IIR filters with high sample rates?

An IFFT of a single non-zero bin (or conjugate mirrored pair) produces a single pure sinusoid (with no phase or amplitude changes).  An IFFT of just a few adjacent non-zero FFT bins produces something like a slowly amplitude modulated sinewave (a beat note) possibly very slowly modulated in frequency.  Note that none of those IFFT results resemble a sinewave that includes a sudden instantaneous phase shift.  The information about the location of a possible sudden phase change is located in all the other FFT bins that an ultra-narrow FFT filter throws away.  e.g. you would have to de-zero a vast number of FFT bins to get that change to reappear in the IFFT result.  That is why sharp changes have very wide-band spectrums.

This renders your proposed filter (zeroing all but a few FFT bins) nearly useless for determining the location of any sharp phase shift.


Interpreting dataset of a power spectrum

1 second at 10ms "hops" gives 100 hops, which correspond to the rows.

In DSP parlance, the frame (the 25ms of samples) is called window size, and is the sample interval that is fed into the DFT (FFT). These samples are then multiplied by a window function prior to calculating the DFT.

25ms with $f_s = 16 \rm{kHz}$ gives 400 samples per frame. From your data, the DFT is performed with a vector size of 512, so probably the frame is then 0-padded from 400 to 512 samples, then fed to the DFT.

The DFT will then produce its output as a vector of frequency samples of size 512, where each frequency corresponds to multiples of $f_s/512 = 31.25 \rm{Hz}$. So your first column corresponds to frequency 0 Hz, then second column to 31.25 Hz, third to 62.5 Hz, etc.

Only half of the output is shown; since your signal is real, the other half of the DFT output (from 8000 to 16000 Hz) will be symmetrical to the first half. Also, since the sampling frequency is 16 kHz, then you cannot represent signals beyond 8 kHz.

Why 257 columns? Probably because the last column corresponds to 8 kHz (inclusive), so you have 256 frecuency intervals, therefore 257 frequency samples.

The plot looks like log of the magnitude of the DFT. The exact units you can deduce studying your audio signal. Try yourself with Octave / Matlab, using the spectrogram functionality there.

The DFT of 1 frame gives you 1 row of your spectrogram. Step forward 10ms, get another frame, repeat.

P.S. Sometimes spectrograms are easier to visualize with time on the horizontal axis.


Should the input of a Kalman filter always be a signal and its derivative?

A state variable and its derivative are often included as inputs to a Kalman filter, but this is not required. The essence of the Kalman framework is that the system in question has some internal state that you are trying to estimate. You estimate those state variables based on your measurements of that system's observables over time. In many cases, you can't directly measure the state that you're interested in estimating, but if you know a relationship between your measurements and the internal state variables, you can use the Kalman framework for your problem.

There is a good example of this on the Wikipedia page. In that example, 1-dimensional linear motion of an object is considered. The object's state variables consist of its position versus time and its velocity on the one-dimensional line of movement. The example assumes that the only observable is the object's position versus time; its velocity is not observed directly. Therefore, the filter structure "infers" the velocity estimate based on the position measurements and the known relationship between velocity and position (i.e. $\dot{x_k} \approx \frac{(x_k - x_{k-1})}{\Delta t}$ if acceleration is assumed to be slowly-varying).


Remember that phase is usually only defined on the $[0;2\pi[$ interval (or on $[-\pi;+\pi[$). When you do that $\mod 2\pi$ operation on your graphs, you'll notice that there's no significant difference between your two plots.


This answer shows how to create a real-coefficient infinite-impulse-response (IIR) filter, the output of which equals the real part of the output of a given complex-coefficient IIR filter. Also an example is given that shows that the approach can be useful.

Z-transform $\mathcal{Z}\{h[n]\}$ of a time-domain sequence $h[n]$ has the property that:

$$\mathcal{Z}\{h^*[n]\} = H^*(z^*),\quad\text{where }H(z) = \mathcal{Z}\{h[n]\}.\tag{1}$$

Given an IIR filter transfer function:

$$H(z) = H_0\frac{\prod^M_{k=1}\left(1-q_kz^{-1}\right)}{\prod^N_{k=1}\left(1-p_kz^{-1}\right)},\tag{2}$$

with real coefficient $H_0$ and $M$ complex zeros $q_k$ and $N$ complex poles $p_k,$ discarding of the imaginary part of its output $h[n]$ to form real output $g[n]$:

$$g[n] = \operatorname{Re}(h[n]) = \frac{1}{2}\left(h[n]+h^*[n]\right),\tag{3}$$

is equivalent to having a filter with transfer function:

$$G(z) = \frac{1}{2}\big(H(z) + H^*(z^*)\big).\tag{4}$$

Equivalently, with $z = e^{i\omega},$ this new filter has a frequency response $\frac{1}{2}\big(F(\omega) + F^*(-\omega)\big),$ where $F(\omega)$ is the frequency response of the full complex filter. Using further properties of complex conjugation:

$$G(z) = \frac{1}{2}H_0\left(\frac{\prod^M_{k=1}\left(1-q_kz^{-1}\right)}{\prod^N_{k=1}\left(1-p_kz^{-1}\right)} + \frac{\prod^M_{k=1}\left(1-q^*_kz^{-1}\right)}{\prod^N_{k=1}\left(1-p^*_kz^{-1}\right)}\right)\\
= \frac{\frac{1}{2}H_0\left(\prod^M_{k=1}\left(1-q_kz^{-1}\right)\prod^N_{k=1}\left(1-p^*_kz^{-1}\right) + \prod^M_{k=1}\left(1-q^*_kz^{-1}\right)\prod^N_{k=1}\left(1-p_kz^{-1}\right)\right)}{\prod^N_{k=1}\left(1-p_kz^{-1}\right)\prod^N_{k=1}\left(1-p^*_kz^{-1}\right)}\\
= \frac{\mathcal{Z}\bigg\{\operatorname{Re}\Big(\mathcal{Z}^{-1}\left\{H_0\prod^M_{k=1}\left(1-q_kz^{-1}\right)\prod^N_{k=1}\left(1-p^*_kz^{-1}\right)\right\}\Big)\bigg\}}{\prod^N_{k=1}\Big(\left(1-p_kz^{-1}\right)\left(1-p^*_kz^{-1}\right)\Big)}.$$$\tag{5}$

By expanding the products to polynomials of $z^{-1}$, the above gives the procedure to calculate the coefficients of the real filter given the coefficients of the original complex filter, by convolutions of sequences of polynomial coefficients. If the original complex filter with input $x[n]$ and output $y[n]$ is in form:

$$\begin{align}y[n] &= b_0\,x[n] + b_1\,x[n-1] + b_2\,x[n-2] + \ldots\\
         &-a_1\,y[n-1] - a_2\,y[n-1] - \ldots,\end{align}\tag{6}$$

that readily accepts the coefficients $B = [b_0,\,b_1,\,\ldots]$ and $A = [a_0,\,a_1,\,\ldots]$ of the numerator and denominator polynomials $b_0+b_1\,z^{-1}+b_2\,z^{-2}+\ldots$ and $a_0+a_1\,z^{-1}+a_2\,z^{-2}+\ldots,$ respectively, with a hidden coefficient $a_0 = 1$, then the coefficients $B_\text{real}$ and $A_\text{real}$ for the new filter will be given by:

$$\begin{align}B_\text{real} &= \frac{1}{2}\left(B*A^* + B^**A\right) = \operatorname{Re}\left( B*A^*\right)\\
A_\text{real} &= A*A^* = \operatorname{Re}\left(A*A^*\right),\end{align}\tag{7}$$

where $*$ represents convolution (or complex conjugate if superscript). The last $\operatorname{Re}$ is unnecessary, but may be useful in numerical convolution.

Example 1: Minimal example

A minimal example in Octave starts with a complex filter and generates a real filter that matches the real output of the complex filter:

b = [0.5, 0.35 - 0.35*j]; #Coefficients of numerator polynomial
a = [1, - 0.6 - 0.6*j]; #Coefficients of denominator polynomial
length = 1024;
impulse = eye(1, length);
[h_complex, w_complex] = freqz(b, a, length, "whole");
[h_discard_imag, w_discard_imag] = freqz(real(filter(b, a, impulse)), [1], length, "whole");
b_real = real(conv(b, conj(a)))
a_real = real(conv(a, conj(a)))
[h_real, w_real] = freqz(filter(b_real, a_real, impulse), [1], length, "whole");
freqz_plot([w_complex, w_discard_imag, w_real], [h_complex, h_discard_imag, h_real]);

The generated real coefficients of the numerator and denominator polynomials of the new real filter are:

b_real =
   0.50000   0.05000   0.00000

a_real =
   1.00000  -1.20000   0.72000

The frequency response of the new real filter equals that calculated from the real part of the impulse response of the complex filter:

Figure 1. Blue: Frequency response of the complex filter with the imaginary part of its output still present. The frequency response lacks Hermitian symmetry, so no real-coefficient filter can reproduce it. Red: Shared frequency response of the new real-coefficient filter and the complex-coefficient filter with the imaginary part of its output discarded. The shared frequency response has Hermitian symmetry.

Example 2: Sliding discrete Fourier transform

Constructing a real filter from a complex filter can sometimes be useful. One such case is calculation of the real part of a single term of a sliding discrete Fourier transform (sliding DFT) of length $r$ at frequency $\omega = 2\pi m/r,\, m\in N.$ The complex IIR filter is:

$$\begin{align}y[n] &= x[n] - x[n-r]\\
&+e^{-i\omega}y[n-1]\end{align}\tag{8}$$

We get using Eqs. 6 and 7 the coefficients of the numerator and denominator polynomials of the transfer function of the real filter $(A_\text{real},\,B_\text{real})$ from those of the complex filter $(A,\,B).$ Assuming that $r > 1$ (coefficients not given are zero):

$$\begin{array}{rcllll}
k&=&0&1&r&r+1\\
\hline
B[k]&=&1&&-1\\
A[k]&=&1&-e^{-i\omega}\\
B_\text{real}[k]&=&\operatorname{Re}(1)&\operatorname{Re}(-e^{i\omega})&\operatorname{Re}(-1)&\operatorname{Re}(e^{i\omega})\\
&=&1&-\cos(\omega)&-1&\cos(\omega)\\
A_\text{real}[k]&=&1&-e^{i\omega}-e^{-i\omega}&-e^{-i\omega}(-e^{i\omega})\\
&=&1&-2\cos(\omega)&1\end{array}$$

The real filter is thus:

$$\begin{align}y[n] &= x[n] -\cos(\omega)x[n-1] -x[n-r] -\cos(\omega)x[n-r-1]\\
&+ 2\cos(\omega)y[n-1] - y[n-2]\end{align}$$

The recursive part can be recognized as the Goertzel filter.

In Octave:

m = 2;
r = 20;
w = 2*pi*m/r;
b = zeros(1, r+2);
b(1) = 1; b(2) = -cos(w); b(r+1) = -1; b(r+2) = cos(w);
a = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
plot([0:length-1], filter(b, a, impulse), "o");

Figure 2. Impulse response of a real IIR filter that calculates the real part of a single frequency bin of sliding DFT.

The impulse response can be zero phase aligned to the opposite end by multiplying the complex transfer function by $e^{-i\omega}:$

m = 2;
r = 20;
w = 2*pi*m/r;
b = zeros(1, r+2);
b(1) = cos(w); b(2) = -1; b(r+1) = -cos(w); b(r+2) = 1;
a = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
plot([0:length-1], filter(b, a, impulse), "o");

Figure 3. Impulse response of a real IIR filter that calculates the real part of a single frequency bin of sliding DFT, with zero phase aligned to the opposite end than in Fig. 2. This is a circular shift of the impulse response one sampling period to the left, and can sometimes be used to remove one sample of system latency.

In case of $m=0,$ the numerator and denominator polynomials will have an unnecessary common factor, but for other values of $m$ the approach works fine.


Hilbert transform of sinusoid -- apparent contradiction

The error lies in the assumption that if $g(t)$ is the Hilbert transform of $f(t)$, then the Hilbert transform of $f(-t)$ must be $g(-t)$. This is not the case.

Let $f^-(t)=f(-t)$. Then we have

$$g(t)=\mathcal{H}\{f\}(t)=\frac{1}{\pi}\text{p.v.}\int_{-\infty}^{\infty}\frac{f(\tau)}{t-\tau}d\tau\tag{1}$$

and

$$\begin{align}\mathcal{H}\{f^-\}(t)&=\frac{1}{\pi}\text{p.v.}\int_{-\infty}^{\infty}\frac{f(-\tau)}{t-\tau}d\tau\\&=\frac{1}{\pi}\text{p.v.}\int_{-\infty}^{\infty}\frac{f(\tau)}{t+\tau}d\tau\\&=-\mathcal{H}\{f\}(-t)\\&=-g(-t)\end{align}\tag{2}$$

So if $f(t)=\cos(\omega t)$ and $g(t)=\mathcal{H}\{f\}(t)=\sin(\omega t)$, then the Hilbert transform of $f^-(t)=\cos(-\omega t)$ equals $-g(-t)=-\sin(-\omega t)=\sin(\omega t)$. Happiness again.

EDIT:

In response to your comment about substituting for $\omega$ instead of reversing time, note that

$$\mathcal{H}\{\cos(\omega_0 t)\}=\sin(\omega_0 t)$$

is only valid for $\omega_0>0$. This can most easily be seen in the frequency domain, where the Hilbert transform corresponds to a multiplication with $-j\;\text{sign}(\omega)$. The Fourier transform of $f(t)=\cos(\omega_0t)$ is given by

$$F(\omega)=\pi\left[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)\right]$$

If $\omega_0>0$, multiplication by $-j\;\text{sign}(\omega)$ gives the Fourier transform of $\sin(\omega_0t)$. However, if $\omega_0<0$, we obtain the Fourier transform of $-\sin(\omega_0t)$, because then $\delta(\omega+\omega_0)$ appears at positive frequencies instead of $\delta(\omega-\omega_0)$.


There appears to be a contradiction with the recurrence solutiom

From the comments under Tim Wescott's answer, it seems that there is still  some confusion concerning the correct way to solve such problems. Let's consider the non-homogeneous difference equation (DE)
$$y[n]-p y[n-1]=q^n,\qquad p\neq q\tag{1}$$
with initial value $y[-1]=c$. Note that initially we don't restrict the forcing function on the right-hand side of $(1)$ to equal zero for $n<0$.
As correctly shown in the OP, the general solution to the homogeneous DE has the form
$$y_h[n]=Ap^n\tag{2}$$
which is also called the natural response,
and a particular solution (forced response) of $(1)$ is given by
$$y_p[n]=Bq^n\tag{3}$$
The constant $B$ is determined by substituting $(3)$ into $(1)$:
$$Bq^n-pBq^{n-1}=q^n\tag{4}$$
from which we obtain $B=\frac{q}{q-p}$. Consequently, the general solution of the non-homogeneous DE $(1)$ is
$$y[n]=y_h[n]+y_p[n]=Ap^n+\frac{q}{q-p}q^n\tag{5}$$
Note that the solution $(5)$ is valid for all values of $n$.
The value of $A$ is determined by the initial condition $y[-1]=c$. Substituting $n=-1$ into $(5)$ and setting it equal to $c$ results in
$$A=cp-\frac{p}{q-p}\tag{6}$$
The solution to the initial value problem is hence given by
$$y[n]=cp^{n+1}+\frac{1}{q-p}\left[q^{n+1}-p^{n+1}\right]\tag{7}$$
We now only consider values $n\ge 0$, simply because that's the range we're interested in.
Note that the term $cp^{n+1}$ in $(7)$ is the zero-input response (ZIR), and the remaining term is the zero-state response (ZSR).
From an engineering point of view, composing the solution into a ZSR and a ZIR appears more natural and more useful than the classical decomposition into a homogeneous and a particular solution, because the decomposition of the solution into ZSR and ZIR separates the contributions arising from an external input and from internal conditions, respectively.


Is the Laplacian Filter an High Pass Filter (HPF)?

There are many approximations for the Laplacian Filter (See The Hypermedia Image Processing Reference - Laplacian/Laplacian of Gaussian):

Indeed this is an High Pass Filter. Namely it will remove low frequencies (Specifically it will remove the DC Value, namely the output image will have mean value of 0).


First, learn about digital differentiators and how to implement them with software. Then apply your accelerometer output signal to a digital differentiator. Next, use software to detect a large positive-valued sample followed by a large negative-valued sample in the differentiator's output sequence. That sample will correspond to your "first selected area". Finally, use software to detect a large negative-valued sample followed by a large positive-valued sample in the differentiator's output sequence. That sample will correspond to your "second selected area".


Map the Amplitude Envelope from [dbmin,dbmax] to [-inf,dbmax]

May be you can try the following:

 dbnew = dbmax * (1 + log((db-dbmin)/(dbmax-dbmin)))

dbnew = -inf when db = dbmin and dbnew = dbmax when db = dbmax. Please let me know if it works.


You are missing one important piece of information in your question and that is the prototype definition for the InPlaceFFT() function. The most important thing to be aware of is that the FFT returns complex valued numbers. It can be that the function expects the data to be formatted such that real and imaginary values are interleaved or it could be that the function expects the real parts to be placed first in the buffer and the imaginary part to be placed in the end of the buffer. It could also be that the function returns data bit reversed (see section "Data reordering, bit reversal, and in-place algorithms" here http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm).

Another potential issue is the length of the buffer you pass as argument, your temp buffer. You might have to pass a buffer with a length that is at least 1024*sizeof(float) or 2048*sizeof(float), if you are sure that the InPlaceFFT() function expects a float* as input. If the InPlaceFFT() function is constrained to work on real input, a buffer size of 1024 might be sufficient. But you should take a better look at the documentation for how the InPlaceFFT() works in order to find out how large the size of you temp buffer should be and how the output is laid out in the temp buffer after it returns. Another thing that looks strange is that you don't need to pass the FFT size to InPlaceFFT().

If you are sure that the sample rate is always 44.1kHz it is fine to hardcode 44.1kHz into the code, otherwise you might want to make your code more flexible by loading the sample rate into a variable and use that in your code.

A 1024 point FFT returns 1024 values (complex). You only need the first 513 output values (these values are typically called bins), the rest are in principle redundant (although the can be handy to have). If you then want to compute the magnitude spectrum you need to compute the absolute value of the complex numbers (so your abs() function should work on complex numbers) and you are in principle done. The only thing left to do is to correct the values from index 1 to index 511 by multiplying them with 2. The values at index 0 and index 512 are real and need no correction. If you are going to plot the magnitude spectrum in a GUI you probably want to plot the magnitudes in dB. So you need to compute the log (base 10) of the magnitude values and multiply by 20.

It would be useful if you could update your post with the documentation for InPlaceFFT().


What you are describing refers to Pitch Shifting.

A naive approach to pitch shifting is based on the shifting property of the Fourier Transform (and more formally here). According to the shifting property, multiplying a signal with a sinusoid in the time domain results in shifting in the frequency domain.

But, convincing pitch shifting takes a little bit more effort than that because it does not only matter to shift the frequencies of part of the spectrum to a new range, it is also necessary to preserve their relative phase for the sound to be recognisable. If you are trying to pitch shift musical instruments (or vocals), then it gets a little bit more challenging, because of the way pitch is working. In that case, you not only have to shift the components in the frequency domain but also preserve their intervals. In other words, it is not only enough to shift those formants that you show in your illustration, you also need to make sure that the new frequencies they end up into, have the same "distances" between them as the ones they are being pitch shifted from. This is not as straightforward as it...sounds.

Pitch shifting in this way implies a linear pitch scale. But the pitch scale is non-linear and going up in pitch (higher octaves) does not result in equally proportional frequencies.

Because of these reasons, more complex pitch shifting algorithms contain a "processing stage" between analysing and resynthesizing the pitch shifted version of the signal in their input. This processing stage attempts to take into account these relationships between the components of the spectrum according to some approach. For example, directly working out these relationships to the new range as the phase vocoder does, or attempting to fit a model to the spectrum and use it to reproduce it at a different range as spectral modelling does.

Hope this helps.


Other techniques like overlap save/overlap add

The traditional overlap-add technique cannot be used as is in musical applications (reverberation) where latency is critical. It is a common practice to break down the impulse response into non-uniform chunks, and to run in parallel different convolvers for each segment of the impulse response - some of them being naive FIR implementations (for the short blocks at the "head" of the response), and some of them being high-latency overlap-add FFT convolvers.

See this paper.


What is the difference between Goertzel algorithm and cross correlation in time

The Goertzel algorithm gets away without calculating a complex product with $e^{j\pi lk}$ n times; thus drastically reducing the computational complexity.


What is the overall capacity of two cascaded BSC channel?

No, this is not correct. Consider the chain of two BSCs with error probabilities $p_1$, $p_2$ as a single BSC with unknown error probability $p$.

Now, we know that in overall no error occurs, when:

neither BSC1 nor BSC2 create an error: $(1-p_1)(1-p_2)$
both channels make an error: $p_1p_2$.

Hence the overall probability of an error is

$$ p= 1-[(1-p_1)(1-p_2)+p_1p_2]=p_1+p_2-2p_1p_2 $$

which is also known as the binary convolution. Now, the capacity of the chain-BSC is given by
$$C=1-H(p)$$
where $H(p)$ is the binary entropy function.


Is the ADC sample-and-hold fast enough to even capture any envelope peaks or even half cycles?  If not, all bets are off.

If your ADC does have a fast enough capture time, then randomizing your sample times might be a good bet at capturing a few near extrema.  If you know, a priori, the approximate shape of your waveform, then Monte-Carlo sampling might allow one to estimate the waveform amplitude from the statistical distribution of enough sufficiently time-randomized samples.  Assuming the envelope changes slowly enough.


Getting 16 bit image statistics in Photoshop

A lot of tools do this by default (even matlab's image processing tools for the most part). I'd load the image up in matlab and use the hist, std, mean, etc function by treating the image as a vector. Tools like EnvI also handle 16 bit images properly.

The matlab/octave code would be roughly as follows (assuming 16 bit greyscale png):

img=double(imread('image.png'));
mean(img(:))
std(img(:))
hist(img(:),number of bins)


Why use transfer functions than differential equations?

Why we use Transfer Functions, when we can get a system's output by just solving it's differential equation?

Because differential equations are unwieldy and hard to deal with, and you can't see the behaviour on different frequencies from these, whereas transfer functions just give you the behaviour of an LTI system given an excitation of given property. Much easier to work with!


Beware.  A Gaussian random walk has a non-zero distribution width that increases with time.  Thus the double integral of that random walk will likely very quickly zoom off the edge of your solution space.  A linear filtering of that random walk won't behave any better.  Whether or not mixed with the "real" acceleration signal.


Since the length of your FFT is exactly one period of your sinusoidal input signal (N = Fs/F0) , the signal will appear in bin 1 (the one above DC) of the FFT result (the conjugate mirror is redundant for strictly real input).  To estimate the phase of the sinusoid of frequency F0, referenced to the start of the data window, just compute the atan2() of the signed imaginary and real components from that single complex bin 1.  No need to abs() components or sum or window or estimate the amplitude or examine any other bin.  Unless you want to sanity check the magnitude of bin 1 to make sure it isn't zero or just numerical noise.


What is the rule for manipulating the boundaries of a summation?

Summations are not like integrals. When you calculate an integral, swapping the boundaries does change the sign of the result. However, when talking about summations:

$$ \sum_{n=-\infty}^{-1}a(n)=\sum_{n=-1}^{-\infty}a(n)$$

This holds because it is the same if we perform the sum of $a(-1)+a(-2)+...+a(-k)$ for any $k>0$ in any order. Namely, you can order those terms as you wish as long as all the terms are present. For infinite summations, this is true if the series is absolutely convergent.

In your example, they just made the change of variables $m=-n-1$ to get the summation from $0$ to $+\infty$. It may be confusing because they kept the original name of the index ($n$) instad of replacing it with the new one ($m$), but remember that they are just indexes or dummy variables, so in the end it doesn't matter. You can check that if $m=-n-1=-(n+1)$, then the boundaries change:

$$n=-1 \to m=-(-1+1)=0$$
$$n=-\infty \to m=-(-\infty)=\infty$$

where I used the fact that a finite number is negligible against infinity. So the summation can be written as:

$$\begin{align} -\sum_{n=-\infty}^{-1} \alpha^nz^{-n}&=-\sum_{m=\infty}^{0} \alpha^{-(m+1)}z^{-(-(m+1))}\\
&=-\sum_{m=0}^{\infty} \alpha^{-(m+1)}z^{-(-(m+1))}\\
&=-\sum_{m=0}^{\infty} \alpha^{-(m+1)}z^{m+1}\\
&=-\sum_{m=0}^{\infty} (\alpha^{-1}z)^{m+1}
\end{align}$$

And that's the same result you wrote in your question.


How do I compute the gradient vector of pixels in an image?

So, if you're using Matlab, you can do:

X = your matrix

[gx,gy] = gradient(X); % first order gradient

[gxx,gxy] = gradient(gx); % second order gradient

[gxy,gyy] = gradient(gy); % second order gradient

To find the curvature of features I advise you to look into the eigenvalues and eigenvectors of a Hessian matrix. A hessian matrix is a square matrix of second-order derivatives of a scalar valued function. In this case the scalar field is the intensity, and the second order derivatives are $g_{xx}$, $g_{xy}$, $g_{yy}$.


If your signal is an analytical signal, then envelope(real(X)) would simply be your answer. This would be identical to using the Hilbert transform to extract the Quadrature part of the signal from the Inphase part. If your signal is not analytical, your envelope would simply be:
$$G = \sqrt{I^2 + Q^2} $$


What kind of edge detection use paint.net?

The edges doesn't have single pixel width, so it's probably not canny (no non-maxima supression). Could be a simple high-pass filter, whose pass-band is masked out in the Fourier domain in specific angles.


Your integrator is a perfect accumulator, which is an unstable (marginally stable) system with a pole at $z=1$:

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

This means that the value of the frequency response at $\omega=0$ (corresponding to the pole at $z=1$) is undefined. Consequently, the value of the group delay at $\omega=0$ is also undefined. Note that you can write down the frequency response of an ideal accumulator using a Dirac delta impulse (see this question and its answer), but this doesn't help in finding a value of the frequency response (or the group delay) at $\omega=0$.

For $\omega\neq 0$ (or, more generally, $\omega\neq 2\pi k$), we can easily determine the frequency response and group delay of the system:

$$H(e^{j\omega})=\frac{1}{1-e^{-j\omega}}=e^{j\omega/2}\frac{1}{e^{j\omega/2}-e^{-j\omega/2}}=e^{j\omega/2}\frac{1}{2j\sin(\omega/2)},\quad\omega\neq 2\pi k\tag{2}$$

From $(2)$, we see that the phase $\phi(\omega)$ equals $\omega/2$ plus a constant term which is irrelevant for the group delay. Consequently, the group delay is given by

$$\tau_g(\omega)=-\frac{d\phi(\omega)}{d\omega}=-\frac12,\quad\omega\neq 2\pi k\tag{3}$$

which is what you found using Matlab. However, Matlab doesn't tell you what's going on at $\omega=0$, where the value of the group delay is indeed undefined.


Frame acquisition vs symbol timing synchronization

Auto-correlation will be enough to determine the symbol center of 1 symbol at the end of the preamble, so that is sufficient for acquisition.

Symbol center tracking throughout the frame is still required to keep the receiver symbol clock synchronized to the transmitter's symbol clock.  Your receiver's symbol clock frequency will never be exactly the same as the transmitter's symbol clock frequency in real implementations, therefore continual estimation and correction or tracking is needed.


Understanding time delay estimation for GPS?

For GPS, (simplifying for now by omitting corrections for ionsphere and orbit position and relativistic clock offsets), we determine the "Pseudo-range" to each satellite (SV), which will be the relative delay between all the received satellites we have correlated to relative to our local clock- using correlation as you described (delay each locally generated SV until you get a maximum peak and that delay value is your pseudo-range. Thus, as long as we have a pseudorange for at least 4 satellites, we can solve a matrix involving 4 equations with 4 unknowns: x position, y position, z position and time (which is the time error between our local clock and the clock of the satellite). The satellite's time keeping is highly accurate (using Atomic clocks) and the 4 equation solution allows us to use very inexpensive clocks in our local receivers.

Below is shown a simplified 2D explanation showing how unique x,y position is determined on a plane with 3 satellites and perfect knowledge of time. Note how with 2 satellites there are two solutions. The same thing occurs with 4 satellites and unknown time, it is just that one of the solutions is on earth and the other can be neglected as it would be out in outer space.  In practice we benefit in a lower error solution by averaging over all satellites in view (forming a least mean square solution to an over determined set of equations).

Below is a sample output of what the output of a sliding correlator looks like for a particular SV in an actual captured GPS signal (the signal itself is indistinguishable from noise since the GPS signal prior to correlation is approximately 20 dB below the noise floor). Here we see the result of the spreading code that repeats every 1 ms, and the specific location of the peak in time relative to any local time reference we choose is the "pseudo-range".  These peaks are what you would get out of the Xcorr command (shown is the absolute value) and the horizontal axis where a peak occurs is the delay you are interested in finding. What is important is that all you need is the relative delay between all the satellites in view and your local clock reference (so you can pick any reference point you want as long as all codes are done with the same clock).  What I have seen in implementation, and perhaps this would be more complicated than you need, is to have a code NCO in a code tracking loop that is used to generate the codes for each satellite, and the code tracking is done using three correlators for each SV in an early-prompt-late arrangement each offset by 1/2 a chip. Early minus Late correlator output would then be proportional to code error (once acquired) and drives the Code NCO to track the delay. In such a fashion the code generator is always kept in alignment with the received code (in contrast to the sliding correlator approach I show below), and the relative delay is determined by the relative position of each code.

Below is the received signal prior to correlation (for the same signal used to generate the correlation plot above), with the blue representing the captured waveform after I had down-converted it to baseband and aligned it with the spreading code shown in red (as evidenced by the correlation peak). To note that any resemblance to the code that you think you might be able to see in the blue waveform to the red waveform is purely coincidental. The blue waveform in this view at this scale is visually indistinguishable from the AWGN process that is dominating it.

This is how the transmitted signal looks at that satellite (just two states +1 and -1):

And this is how that same signal looks as captured at our receiver, after down converting to baseband, removing all carrier and code offsets, and sampling at the mid point of each symbol (to match what we know was +1/-1 at the satellite). This matches very well to an AWGN noise process within the sampling space described, with no evidence of the signal without performing the correlation with the known sequence.

To add to the simplified descriptions above, here is a link to a good reference on the complete PVT solution (position, velocity and time) using GPS: http://www.nbmg.unr.edu/staff/pdfs/Blewitt%20Basics%20of%20gps.pdf


Is there any fast algorithm to compute the most similar rectangular regions in two images?

As stated already, you can not do much to speed up a comparison between any two chosen patches. What you need to focus on in order to speed up the process is how to reduce the number of patch pairs you need to compare.

If the computational process is as expensive as I presume it is, in addition to already mentioned image pyramid, I might have another suggestion.

Extract local features and perform feature matching. This is the way to usually recognize similar images, but it works for finding objects in whole images, which means it should also match features in smaller, locally similar patches.

These guys: Sivic, Zisserman: Video Google: A text retrieval approach to object matching in videos had a good idea about spatial consistency. In short, they make sure that the groups of matches have a fairly similar geometrical layout. They use a fairly loose criterion but suggest a way to both loosen and strengthen the criterion. You do have to do a straightforward brute force matching between the feature descriptors, but that still seems far less computationally expensive than directly calculating similarities between patches.

This way, you will hopefully get clusters of matches. Then you could limit the search only to the patches that contain a certain minimal number of matches. In my opinion, these patches would have a substantially higher probability of being the most similar. This is a bit complex add on to make, but I think it could be worth a shot.


Fourier Like Spectral Analysis with Uneven Intervals and Redesigned DFT Matrix

I have derived the DFT for data whihc is sampled in a non uniform manner:
The DFT Matrix for Non Uniform Time Samples Series
Problem Statement
We have a signal $ x \left( t \right) $ defined on the interval $ \left[ {T}_{1}, {T}_{2} \right] $.
Assume we have $ N $ samples of it given by $ \left\{ x \left( {t}_{i} \right) \right\}_{i = 0}^{N - 1} $. The samples time $ {t}_{i} $ is arbitrary and not necessarily uniform.
We're after the DFT of the samples $ \left\{ X \left[ k \right] \right\}_{k = 0}^{K - 1} $ as it was samples in a uniform manner (Implicitly means the samples in Frequency domain will be uniform as well).
Deriving the Connection
In the DFT Transform the connection between time and frequency is given by:
$$ x \left[ n \right] = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N} n } \tag{1}  $$
In $ \eqref{EqnIdft} $ we use $ n $ for modeling the sample index in time. We usually build samples in time as $ x \left[ n \right] = x \left( n {T}_{s} \right) $ where $ {T}_{s} $ is a uniform sampling interval.
Hence we could write:
$$ x \left( n {T}_{s} \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} n {T}_{s}} \tag{2} $$
In $ \eqref{EqnIdft2} $ we added explicit scaling of time. This is a known property of Fourier transform family which scales the domain in order to normalize the transform.
Now, there is nothing which blocks us from using arbitrary time:
$$\begin{align*} \tag{3}
x \left( t \right) & = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} t} && \text{} \\
& = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k {F}_{s}}{N} t} && \text{Since $ {F}_{s} = \frac{1}{{T}_{s}} $}
\end{align*}$$
As can be seen $ \eqref{EqnIdft3} $ makes sense as it goes through each element according to its frequency and sums to give the output at time $ t $. We can go step farther and generalize it for cases we don't have uniform sampling frequency.
The average sampling frequency is given by $ \bar{F}_{s} = \frac{N}{ {T}_{2} - {T}_{1} } $. Let's define $ T = {T}_{2} - {T}_{1} $ and we'll get:
$$ x \left( t \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} } $$
Which is many ways resembles the DTFT Transform equation which does the same in the other direction, transforming uniform discrete samples in time domain to arbitrary frequency (Within a frequency interval) in Frequency Domain:
$$\begin{align*} \tag{4}
X \left( f \right) & = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi f {T}_{s} n } && \text{} \\
& = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi \frac{f}{ {F}_{s} } n } && \text{Since $ {F}_{s} = \frac{1}{{T}_{s}} $}
\end{align*}$$
We see the same scaling, $ \frac{f}{ {F}_{s} } $ which scales the continuous $ f $ relative to the interval of frequencies $ {F}_{s} $ which is equivalent to $ \frac{t}{ T } $ which scales $ t $ relative to the time interval of the continuous signal.
The Transform Matrix
So, given the set of time indices $ {\left\{ {t}_{i} \right\}}_{i = 0}^{N - 1} $ the transformation matrix, from frequency domain to time domain, is given by:
$$ D \in \mathbb{R}^{N \times K}, \;  {D}_{i, k} = {e}^{ j 2 \pi k \frac{ {t}_{i} }{T} } $$
The Model
In vector form the model is:
$$ x = D y $$
Where $ y \in \mathbb{C}^{K} $ is the vector of the frequency coefficients in uniform grid, $ x $ is the samples in time (Non Uniform, Or at least no assumption of uniformity) and $ D $ as defined above.
Since in our model we're after $ y $ the answer is given by:
$$ y = {D}^{\dagger} x $$
Where $ {D}^{\dagger} $ is the Pseudo Inverse Matrix of $ D $.
Implementation & Results
The code is as following:
subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

samplingFrequency = 101; %<! [Hz]
samplingInterval = 1 / samplingFrequency; %<! [Sec]
startTime = 1; %<! [Sec]
endTime = 4; %<! [Sec]
timeInterval = endTime - startTime; %<! [Sec]

numSamples = round(samplingFrequency * timeInterval);
numSamplesTT = round(1.2 * numSamples);

signalFreq = 2; %!< [Hz]

% The uniform time grid
vT      = linspace(startTime, endTime, numSamples + 1);
vT(end) = [];
vT      = vT(:);

% The non uniform time grid - Reconstruction
vTT = endTime * rand(numSamplesTT, 1);
vTT = sort(vTT, 'ascend');

% The non uniform time grid - DFT
vTD = linspace(startTime, endTime, (10 * numSamples) + 1);
vTD(end) = [];
vTD = vTD(sort(randperm(length(vTD), numSamples)));
vTD = vTD(:);

% The uniform frequency grid
vF      = (samplingFrequency / 2) * linspace(-1, 1, numSamples + 1);
vF(end) = [];
vF      = vF(:);

vK = [-floor(numSamples / 2):floor((numSamples - 1) / 2)];
vK = vK(:);

%% Generate Data

vX  = cos(2 * pi * signalFreq * vT);
vFx = fftshift(fft(vX));

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = subplot(1, 2, 1);
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['Reference Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);

hAxes           = subplot(1, 2, 2);
hStemObj = stem(vF, abs(vFx));
set(hStemObj, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['DFT of the Reference Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
    'FontSize', fontSizeTitle);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - Reconstruction

mD = exp(1j * 2 * pi * (vTT / timeInterval) * vK.') / numSamples;

% Reconstruction according to the model
vY = real(mD * vFx);

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTT, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - DFT of the Non Uniformly Sampled Data

vY  = cos(2 * pi * signalFreq * vTD);

mD = exp(1j * 2 * pi * (vTD / timeInterval) * vK.') / numSamples;
vFy = pinv(mD) * vY;

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTD, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
set(hAxes, 'NextPlot', 'add');
hStemObj    = stem(vF, abs([vFx, vFy]));
set(hStemObj, 'LineWidth', lineWidthNormal);
% hLineSeries     = plot(vTT, vY);
% set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['DFT of the Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

Results are:

Summary
In this post we derived how to estimate the Uniform DFT of a Non Uniform Time Series by solving linear system of equations.
The full code is available on my StackExchange Signal Processing Q32137 GitHub Repository.
Remark: Why Do We Need to Apply fftshift() on the DFT of the Signal?
Indeed in the Reconstruction part we use fftshift(). The shallow answer is easy, we also build the vector vK as symmetric around zero.
But there is a deeper reason for that. In the DFT when we use uniform sampling in Frequency Domain and Time Domain Magic happens without us seeing it explicitly.
When we defined the term $ \frac{k}{ N {T}_{s} } n {T}_{s} $ we replaces $ n {T}_{s} $ with $ t $ hence we prevent the term $ {T}_{s} $ to cancel itself. Now setting $ {F}_{s} = \frac{1}{{T}_{s}} $ means that we multiply by $ k $ and we get frequencies which are out of the Nyquist Frequency.
In most cases when we that happens the Modulo property of the exponent comes in and we get the correct negative value of the frequency in the range $ \left[ -\pi, \pi \right] $. Yet when $ t $ is arbitrary we can think that $ {F}_{s} $ is changing per sample which means when we go farther than $ \pi $ the modulo doesn't bring us to the correct answer.
First, as intuition, always think the DFT is defined on the $ \left[ -\pi, \pi \right] $ interval and it is continuous. So as long as you work on this range things works as intended. This intuition can come from the Fourier Series and Discrete Fourier Series (DFS).
Let's try explaining it using a concrete example. Let's examine the exponent term from the derivation:
$$ 2 \pi \frac{k}{N {T}_{s}} n {T}_{S} = 2 \pi \frac{k}{{F}_{S}} \frac{{F}_{s}}{N} n = 2 \pi \frac{k b}{{F}_{s}} n $$
Where $ b $ is the Bin Resolution in the Frequency domain. Now given the signal is:
$$ x \left( t \right) = \cos \left( 2 \pi f t \right) \Rightarrow x \left( n {T}_{s} \right) = \cos \left( 2 \pi f {T}_{s} n \right) \Rightarrow x \left[ n \right] = \cos \left( 2 \pi \frac{f}{ {F}_{s} } n \right) $$
For $ {F}_{s} = 100 $ [Hz] and $ N = 100 $ (Which means $ b = 1 $) we will have delta at $ k = 2 $ and $ k = 98 $. For $ k = 98 $:
$$ 2 \pi \frac{98}{{F}_{s}} n $$
This is clearly above the Nyquist frequency ($ \frac{{F}_{s}}{2} $) and only for $ {F}_{s} = 100 $ its modulo is $ -2 $ which is correct. But in the model above, since we have arbitrary $ t $ one could think we have changing $ {F}_{s} $ which means we don't get the correct value.
This means the actual equation should be:
$$ x \left( t \right) = \frac{1}{N} \sum_{k = \left \lfloor - \frac{N}{2} \right \rfloor }^{ \left \lfloor \frac{N - 1}{2} \right \rfloor } X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} } $$


I believe there's either a mistake in the presentation or the presentation is using a different definition of linear.

For example, the system is linear in $x$ from a system perspective, but it's affine in $x[n]$ (and, therefore not linear) because of the $bx[n-3]$ offset.

On this site, we tend to go with the system definition rather than split hairs about linearity versus affine-ness.


Decomposition of **3D** structuring elements for morphological operations

I guess this depends on the digital distance transform that one is approximating on the 3d grid and there are various local connectivities possible. There is an implementation in ImageJ here.

It would also be good to verify if you are using a non-flat structuring element or a correct 3d structuring element. Read Matlab reference here. In the place of euclidean distance sqrt(x.^2 + y.^2 + z.^2) one could add the Manhattan distance using this. Quick 3d distance here.

For 3d structuring element decomposition one can see. But the 3d decompositions of a convex shapes like sphere into separable 2d lines is non-trivial, certain shapes are easier than others like the cube. One can refer here for efficient algorithms for spheres.


Cyclic convolution between two matrices in matlab

Assuming that cconv2 does not exist, you can do it with the Fourier transform and multiplication:

A = some_MxN_matrix;
B = some_KxP_matrix;
Mfft = max(M,K);
Nfft = max(N,P);
C = ifft2( fft2(A,Mfft,Nfft) .* fft2(B,Mfft,Nfft) );

Hopefully the above helps.


variable length no-refrence filter algorithm

One possible implementation would be to run two (or more) different order filters in parallel, and use a mixer to blend them by ramping the gains up and down depending on your selection criteria.


scipy.signal.decimate vs array operations

You should use the 'full' parameter for the convolution. Then, naturally with the convolution, your peak that you use for decimation occurs at the end of each symbol. At this point, the transmit and receive filter line up and you get the maximum overlap and hence best SNR.
The decimate function applies a low-pass filter before the decimation. This is due to the sampling theorem: When you reduce the sampling rate, you need to remove higher frequencies to prevent aliasing. However, in your case, you have already convolved with your filter, so the simple version of taking only every n-th sample is perfectly correct.


Errors when converting from digital filter object to transfer function

Generally you should never use the Transfer Function syntax. Look how your poles are closely placed next to the unit circle. With high order Butterworth filters you will see the effect of round-off numerical errors. What is the high order? Well 5 is already high and your filter has 30.

Instead you should use SOS or Zero-Pole-Gain representation. For example try to visualise the same filter but with SOS matrix:

fvtool(d.Coefficients) % or: zplane(sos2zp(d.Coefficients))

You will see that there are no visible numerical errors.


What is the definition of small/large signal bandwidth?

In communications, whether a signal is narrowband (has a small bandwidth) or wideband (has a large bandwidth) is relative to the channel's impulse response.

If the channel behaves as a flat filter (constant gain) over the signal's bandwidth, then the signal is narrowband. If the channel introduces significant distortion (that is, if different signal frequencies see different gain) then it is said that the signal is wideband.


Why are vector quantizers not used in video codecs?

Theoretically, VQ is always better than scalar quantization (SQ). However, as you already mentioned, the complexity of VQ is higher than the complexity of SQ. It increases exponentially together with the dimensionality of vectors. Therefore, the question is how to design a practical coding system which can provide a best rate-distortion performance under a given complexity constraint.

A designing of a codebook is another very important issue in VQ. The main practical concerns about using a codebook in VQ are:

How to generate a universal codebook for a large number of images.
How to scale the codebook for various rate-distortion requirements.
How much storage space is reasonable for storing the codebook in memory.

All this practical constraints make it difficult for VQ to force through standardization process.


Understanding Parseval's Theorem with Discrete Wavelet Transform

Parseval's identity and Plancherel's theorem finally boil down to orthogonality. When one decomposes a data (with samples), via a scalar product, onto an  orthogonal sequence (yielding coefficients), there exists a certain preservation (equality, up to a proportionality factor) of energy between samples and coefficients. There are some technical conditions, and in certain cases, one only get inequalities (cf. Bessel's inequality) or frame bounds.
The equation for the discrete wavelet transform (DWT) might be incomplete, with respect to the indices. I for instance think that in the second term of the RHS, the scaling factor should be $N_j$, not $N_J$ (and this depends a bit on how discrete wavelets are implemented). Basically, an orthogonal wavelet transform projects data onto basis elements gathered in groups called subbands. Each wavelet subband comes from $N_j$ vectors, with an additional $N_J$ vectors for the approximation. And normally, the total number of vectors  should be (about, honestly, this depends on signal extension) the number of samples $N$, in other words: $N=N_J +\sum_{j=1}^JN_j$.


Image processing (number recognition) woes

You may want to check PyTesser. Also, this article is worth reading.


Perhaps the easiest way to understand it is by seeing what happens as we plot graphs of an arbitrary sequence $x[n]$:

As it goes through the system and take one every second sample (i.e. $x[0]$, $x[2]$, $x[4]$, $x[6]$, $x[8]$ and $x[10]$ in the graph above), you get $y[n] = x[nM]$ which can be depicted as:

Now when you shift the input you get $x'[n] = x[n-1]$:

If we again take one out of every second sample we would then get $y'[n] = x'[nM]$:

We then note that this corresponds to the sequence $x[1]$, $x[3]$, $x[5]$, $x[7]$ and $x[9]$ from the original graph. For this simple case where $M=2$ you can see that the result is
$$
\begin{align*}
  y'[n] &= x'[nM] \\
        &= x[2n-1] \\
        & =x[nM-1] &\mbox{given $M=2$}
\end{align*}
$$

A similar process could be followed to arrive at the general result whereby the output of that system given an input of $x[n-1]$ would be $x[nM-1]$.


Finding where signal is showing some (unknown) periodicity

If you're gonna use ASDF to measure periodicity and you want to window the data with something other than rectangular window, do it after subtracting and squaring:

$$ Q_x[k, n_0] \triangleq \frac{2}{N} \sum\limits_{n=0}^{N-1} \left(x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor] \ - \ x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor + k] \right)^2 w\left(\tfrac{n}{N}\right) $$

where

$\left\lfloor \cdot \right\rfloor$ is the floor() function and, if $k$ is even then $ \left\lfloor \frac{k}{2}\right\rfloor = \left\lfloor \frac{k+1}{2}\right\rfloor = \frac{k}{2} $

and $w\left(\tfrac{n}{N}\right)$ is a window function of non-zero width $N$ samples centered at $n=\tfrac{N}{2}$.  A Hann window would be

$$ w(u) \triangleq \begin{cases}
\tfrac12 - \tfrac12 \cos(2\pi u) \qquad & 0 \le u < 1 \\
\\
0 & \text{otherwise} \\
\end{cases}$$

To make this ASDF into "autocorrelation" (in the neighborhood of sample $x[n_0]$) defined from the ASDF:

$$ R_x[k,n_0] = R_x[0,n_0] - \tfrac12 Q_x[k, n_0]  $$

where

$$ R_x[0, n_0] \triangleq \frac{2}{N} \sum\limits_{n=0}^{N-1} \Big(x[n+n_0-\left\lfloor \tfrac{N}{2}\right\rfloor]\Big)^2 w\left(\tfrac{n}{N}\right) $$

Since $Q_x[0, n_0] = 0$ and $Q_x[k, n_0] \ge 0$ for all lags $k$, that means that $ R_x[k, n_0] \le R_x[0, n_0] $ for all lags $k$.

Suppose for a minute that $x[n]$ is periodic with period $P$ (and $P$ happens to be an integer), then

$$ x[n+P] = x[n] \quad \forall n $$

and  $Q_x[mP, n_0] = 0$ and $R_x[mP, n_0] = R_x[0, n_0] \ge R_x[k, n_0]$ for any integer number of periods ($m$ is an integer).  So you get a peak at $k=0$ and at $k$ equal to any other multiple of $P$ if $x[n]$ is periodic.  If $x[n]$ is not perfectly periodic, what we might expect is the biggest peak at $k=0$, another peak (but slightly smaller) at $k=P$ (the period we are looking for) and progressively smaller peaks for larger multiples of $P$.

The measure of periodicity in the neighborhood of $n_0$ would be

$$\frac{R_x[P, n_0]}{R_x[0, n_0]}$$

Hmmm, I just thought of something.  Because of my audio/music-centric POV, I have been assuming no DC bias (i.e. $x[n]$ should maybe come out of a DC blocking high-pass filter).  But I don't think that is the case with your application.  So maybe just stick with ASDF, but somewhere you will need a threshold for how low $Q_x[k, n_0]$ can go (for $k \ne 0$) for you to consider that $x[n]$ is periodic in the neighborhood of $n_0$.


Is jointly wss (wide sense stationary) a transitive relation?

The simplest counter example I can think of (and it is certainly pathological) is $y(t)=0$ for any $x,z$ which are not jointly WSS.


Why the spectral coherence is unity for all frequencies between single-frequency time series and itself

According to the documentation, the coherence between two signals $x(t)$ and $y(t)$ is defined as $C_{xy}(w) = \frac{|P_{xy}(w)|^2}{P_x(w) P_y(w)}$, where $P_x(w)$ and $P_y(w)$ are power spectral density estimates of $x$ and $y$ and $P_{xy}(w)$ is an estimate of the cross-spectral density.

In your first experiment, you put in two identical signals $x(t) = y(t) = \cos(2\pi w_1 t)$. Since $x(t)=y(t)$, what you should expect is that $P_x (w)= P_y(w) = P_{xy}(w)$ for all $w$. Therefore, $C_{xy} = \frac{|P_{xy}(w)|^2}{P_x(w) P_y(w)}= 1$ for all $w$.
You might want to argue that the signal(s) only contains a single frequency. In this case we would expect that $P_x(w) =P_y(w) = P_{xy}(w) =0$ for $w \neq w_1$. However:

This would mean that $C_{xy}(w) = \frac{0}{0} = $undefined for $w\neq w_1$.
This wouldn't happen numerically, as due to rounding errors you'd have some epsilons left that could still make this fraction equal to one
You're not dividing power densities but (Welch) estimates of the power densities of a finite window of the signal, hence $P(w)$ will not be exactly zero anywhere. It will contain estimation artifacts and the truncation effects. Since $x(t)$ and $y(t)$ are exactly equal, the artifacts are exactly equal and the coherence is one.

When you add noise to $x(t)$ and $y(t)$ in your second experiment you're adding "signal content" to the frequencies outside $w=w_1$. Since the noise contributions are drawn independently, they are uncorrelated which means that $P_{xy}(w)$ will indeed be zero there, while $P_x(w)$ as well as $P_y(w)$ will not be. Hence the coherence does go to zero.


MATLAB : Proper estimation of weights and how to calculate MSE for QPSK signal for Constant Modulus Algorithm

Your code reveals many misconception about what the CMA is supposed to achieve:

your step size mu is much too small; note, however, that the optimal step size can only be found through experiment.
the variabe noisedB appears to be the desired SNR of the received signal. An SNR of $0\,\text{dB}$ as specified by you is very poor (the noise is as strong as the signal), and in such noisy conditions the CMA cannot perform properly.
you have a channel with $3$ taps, and your equalizer length is also only $3$ taps; this is much too short to achieve any reasonable equalization.
you define the MSE as the error between the channel impulse response and the equalizer response. This is wrong, because the equalizer is not supposed to converge to the channel impulse response! The equalizer must equalize the channel, i.e. the concatenation of the channel and the equalizer should give a pure delay.

There are a few other small errors in your code, but instead of correcting them all I show you a very simple and short code example. The figure below the code shows the error curve vs symbol index, and the total impulse response of the channel and the equalizer. As you can see, it approximates a delta impulse (i.e. no distortion). The SNR is $30\,\text{dB}$ for this example. The error becomes much smaller if there is no noise. For worse SNRs, the algorithm becomes unstable.

N = 20000;          % # symbols
h = [1,.45,-.2];    % channel impulse response
h = h/norm(h);
Le = 20;            % equalizer length
mu = .001;          % step size
snr = 30;           % snr in dB

s0 = round( rand(N,1) )*2 - 1;     % BPSK signal
s = filter(h,1,s0);                % filtered signal

% add Gaussian noise at desired snr
n = randn(N,1);
vs = var(s);
vn = vs*10^(-snr/10);
n = sqrt(vn)*n;
r = s + n;          % received signal

e = zeros(N,1);     % error
w = zeros(Le,1);    % equalizer coefficients
w(Le)=1;            % actual filter taps are flipud(w)!

yd = zeros(N,1);

for i = 1:N-Le,
    x = r(i:Le+i-1);
    y = w'*x;
    yd(i)=y;
    e(i) = y^2 - 1;
    w = w - mu * e(i) * y * x;
end

np = 100;           % # sybmols to plot (last np will be plotted); np < N!

subplot(3,1,1), plot(e.*e), title('error')
subplot(3,1,2), stem(conv(flipud(w),h)), title('equalized channel impulse response')
subplot(3,1,3), plot(1:np,s0(N-np+1:N),1:np,yd(N-np+1-Le+1:N-Le+1))
    title('transmitted and equalized signal'), legend('transmitted','equalized'), axis([0,np,-1.5,2])


Compensate delay of non-symmetric (non-linear phase) filter

If your filters have a low pass characteristic (which at least the one in your example appears to have), then a possibly useful measure for the delay is the center of gravity of the impulse response:

$$\hat{\tau}=\frac{\displaystyle\sum_nnh[n]}{\displaystyle\sum_nh[n]}\tag{1}$$

For more general impulse responses, other than low pass, you can modify $(1)$ by either taking magnitudes or squares of $h[n]$. Note that $(1)$ is the filter's group delay evaluated at DC.


Assuming that you have the correct formula for the BER (which is probably based on the assumption that a Gray code is used), then - as mentioned by MBaz in a comment - you would have to invert the $Q$ function, which is a problem for which no analytic solution exists. You can use numerical methods implemented in Matlab and probably in other software packages, but you can also just estimate the required value of $E_b/N_0$ from a graph of the bit error rate. I would assume that the latter is what you're supposed to do if that's a homework problem.

Simple approximations (actually bounds) of the $Q$ functions exist, but inverting them will very likely result in very inaccurate estimates of $E_b/N_0$.


One natural way to do this is to use the the total $2N$ samples (two chunks of $N$ samples each) to estimate the parameters. Amplitude, frequency, and phase estimation are all textbook equations that you can look up so the problem here is the fact that they are all mixed together. If you are assuming that you don't know how many sinusoids are mixed, then you also have to detect how many are present. The basic pseudo code is:
Detect how many sinusoids there are
For each sinusoid:
    Estimate amplitude
    Estimate frequency
    Estimate phase
end
Reconstruct signal at n = N + 10000

Detect how many sinusoids there are
You can do this by detecting how many peaks there are in the FFT of the signal. A sinusoid will present a peak at its frequency. At a peak, the slope of the curve will change from positive to negative, so the slope has a zero crossing at the peak. For a discrete signal, there are different ways to compute the slope/derivative but I would start with the simple before moving to more sophisticated methods. A simple way to compute the derivative is just the difference, $x_d[n]=x[n]-x[n-1]$, and a zero crossing happens if $x_d[n-1] > 0 > x_d[n]$. This will give you the number of zero crossings and where they occur.
Edit
The model that was presented in the question is a noise free model, and the detection described relies on that. You can still use this zero crossing method for a noisy signal but need to apply a threshold because the noise will cause many small peaks to occur. The threshold is needed to keep only the ones that have a significant enough amount of energy. Threshold should be picked to satisfy some probability of false alarm metric that you desire.
Estimate ...
As I mentioned above, estimating the amplitude, frequency, and phase of a sinusoid can be found through the FFT peak finding method (and we already found the locations of the peaks) so we are all set.
Reconstruct signal
Now that we have all the estimates, you can calculate your equation using your estimates.
$$ \hat{s}(n+10000) = \sum_i \hat{A}_i \text{sin}\bigg( \frac{2\pi \hat{f}_i (n+10000)}{f_s} + \hat{\theta}_i \bigg) $$


How do I find the frequency of a particular bin in a cosine transform, and two other questions

DCTs (they exist under different flavors) map data along some cosine-related frequency axis, from $0$ to... the maximum frequency of your data, which we don't know, since we don't have its sampling frequency. So the realtive frequency index is implicit, between $0$ and $255/256$ in relative frequency. $2.669825 $, mostly probably, is the DC component, here the average of your signal, divided by $\sqrt{N}$, where $N$ is your number of samples ($256$). If I try to reconstruct your data with a type-II DCT, here is what I get:

The so-called "Frequency index" should be taken with caution, as it really depends on the DCT flavor, as commented be RBJ.

I'm not sure the signal is exactly yours, as I did a guess on the DCT type. To filter above $10$ Hz, you can see to zero the amplitude of all frequency indices above some ratio related to your sampling frequency. DCTs are discrete transforms, acting on signals with "unit" sampling period. So without data sampling information, you cannot get Hertz from the DCT spectrum alone, unless you know spectral characteristics on your signal that can be detected in the spectrum.

Negative values simply mean that the cosine contribution at that frequency is negated.

But if you want to filter data, using a discrete Fourier transform is probably more appropriate.


Removing periodic noise from image using Fourier transform

Just eyeballing the first image tells me that you have a very simple band pattern. They're in fact mostly horizontal (in the spatial domain). You see a clear cross in the Fourier domain. That's because most of your spatial features in that hallway are oriented horizontally and vertically. Your horizontal band similarly show up as a pair of dots above and below the origin. The bands are slightly rotated to the right, and so will the pair of dots.

You'll need to erase a very small area around these dots. But how to locate them? The easiest solution may be to take your recorded image, open it in Paint or a similar program, and paint pure black and white bars over the actual bands. Now take the FFT of this bar pattern, and you get the exact location of the two points you're looking for. It in general won't be exactly a single FFT bin, so you could fit a gaussian curve to the neighborhood of those two points, and use the opposite 1-gauss(x-xc,y-yc) as the frequency domain filter.


Optical realization of the Radon Transform

I actually ran into this same paper and had the same problem.

The paper linked below does a good job addressing the issue. Starting on page 13 the discuss the issue with this optical element.

https://mstamenk.github.io/assets/files/OpticalRadon.pdf


Using DFT Rotation for Frequency Down-Mixing

Is it possible to perform Frequency Down-Mixing ... by simply rotating the DFT sequence?

Yes.

If so what are the Advantages and Disadvantages of this method?

As compared to the standard method of just multiplying with a cosine or complex exponential, there are no advantages (hence, no one really does it). There are lot of disadvantages: it requires framing including overlap management and potentially windowing, it's computational expensive and calculates a lot of stuff that you actually don't need.


The Gibbs phenomenon occurs wherever a discontinuous function is approximated by a truncated (Fourier) series. For the design of systems, the Gibbs phenomenon occurs when approximating ideal brick wall filter responses (like low-pass, band-pass, etc.), but also for all-pass systems such as the Hilbert transformer because of its discontinuity at $\omega=0$:

$$H(e^{j\omega})=-j\textrm{ sign}(\omega),\quad |\omega|<\pi\tag{1}$$

The impulse response of the ideal discrete-time Hilbert transformer is

$$h(n)=\left\{\begin{array}{lc}\frac{2\sin^2(\pi n/2)}{\pi n},&n\neq 0\\
0,& n=0\end{array}\right.\tag{2}$$

A naive design approach would be to simply truncate the impulse response given by (2) to obtain an FIR approximation of a Hilbert transformer. The magnitude of the resulting frequency response for an FIR filter with 23 taps is shown by the blue curve in the following figure:

The point '1' on the x-axis is the Nyquist frequency. One can clearly see the Gibbs phenomenon with large ripples near the discontinuity. Note that with an even length filter the approximation would be better around the Nyquist frequency (because such a filter doesn't have a zero at Nyquist), but the Gibbs phenomenon around 0 would remain. The green curve in the figure is the magnitude of the frequency response when the truncated impulse response is windowed by a Hamming window. Clearly, the effect of the Gibbs phenomenon is greatly reduced, but this comes at the expense of a wider transition band.

Using more sophisticated methods for designing FIR Hilbert transformers, such as the Parks McClellan algorithm, weighted least squares, or constrained least squares, the Gibbs phenomenon can be avoided. Of course, there will always remain an approximation error, but it can be made sufficiently small by choosing an appropriate filter length. Note that this is not the case when simply truncating the ideal impulse response. In this case, the maximum ripple size near the discontinuity remains the same, regardless of the filter length. Also note that with FIR filters the phase approximation error can be made exactly zero by choosing an anti-symmetric impulse response. The only error is the error in the magnitude of the response, as shown in above figure. This is different from IIR approximations, which can achieve an exact magnitude of 1, but there will always be a phase approximation error. You can find some basics about the design of Hilbert transformers in the textbook Discrete-time Signal Processing by Oppenheim and Schafer.


Equivalence of difference equation with system of difference equations

HINT:

Just express $y[n]$ and $x[n]$ in terms of $v[n-k]$, $k=0,1,\ldots$, and show the equality.


It's the key point of array signal processing, I suppose. Say $x$ is the input vector of $[N,1]$ dimension collected from $N$ array sensors. $x(k)$ is its realization at the $k$ moment of time. By its definition covariance matrix (sometimes it's called autocorrelation matrix):

$R = E[x\cdot x^H]$ ,

where $E[]$ is expectation operator and $x^H$ is Hermitian conjugate. For the ergodic process

$R = \lim_{M\to\infty} 1/M \cdot \displaystyle\sum_{k=0}^M x(k) \cdot [x(k)]^H$ .

But in practice we can estimate $R$ with necessary precision with the snapshot of finite length. In the array processing theory snapshot is a group of vectors $x(k)$. It's the basic data block for array processing algorithms.

$R = 1/K \cdot \displaystyle\sum_{k=0}^K x(k) \cdot [x(k)]^H$ ,

where $K$ is the number of spatial vectors or snapshot size. The interesting question is what is optimal value for $K$. In models I've done, $K$ varies from 64 to 256. The estimation precision is enough for MVDR beamformer or Capon spectral estimation. There is interesting trick for adaptive beamformer design (if you have further interest). You can reduce $K$ (and accelerate adaptation process) if you use so-called diagonal loading, look:

$R = R + \sigma^2 \cdot I$,

where $\sigma$ is some variable depending on SNR (maybe someone define it more precisely?) and $I$ is identity matrix of the size $[N, N]$, same size as $R$ actually. But if you want to do spectral estimation procedure by your array (DoA estimation), you shouldn't perform diagonal loading because it will rise noise floor of the estimation and some weak signals of interest will be lost.

Covariance matrix is the second order statistic of the random process which is measured at the array sensors. It contains information about the sources in space (number, strength, direction) and can be used for sources detection and separation. Actually the number of independent spatial signals at the input of array and the rank of $R$ is the same. Singular values decomposition (SVD) of $R$ gives the information about signal subspace which is necessary for subspace-based DoA estimation techniques like MUSIC or ESPRIT. It is worth to say algorithms based upon covariance matrix inversion/decomposition suffers from source correlation. If the sources is highly correlated (have the same waveforms or too close direction) they can't be separated. And also the system performance degrades in highly correlated scenario.

Very good reference about the topic will be

  Harry L. Van Trees - Detection, Estimation, and Modulation Theory, Optimum Array Processing

In the book more general form of covariance matrix is discussed. It's shown that it can be defined either in time domain or in frequency domain. And furthermore frequency domain interpretation is more common since it allows us to solve the problem of wideband beamformer (or DoA estimation) very easy.

Hope this helps.


Why is circular convolution used in DSP? Why not linear convolution?

Given a discrete time LTI system with impulse response $h[n]$, one can compute its response to any input $x[n]$ by a convolution sum: $$y[n] = x[n] \star h[n] =  \sum_{k=-\infty}^{\infty} {h[k]x[n-k]} \tag{1}$$

Without anything further stated, above definition is for the linear convolution (aperiodic convolution) between $h[n]$ and $x[n]$, which are aperiodic discrete-time sequences of possibly infinite length, unless otherwise stated.  This is different from a circular convolution which is between two periodic sequences of period $N$, and computed over a single period.

You can compute a linear convolution in time domain by Eq.1, or in frequency domain using the following DTFT (discrete-time Fourier transform) property:
$$y[n] = x[n]\star h[n] \implies  Y(e^{j\omega}) = X(e^{j\omega}) H(e^{j\omega}) \tag{2}$$

DTFT is naturally related with the linear convolution, as it deals with theoretically existing aperiodic sequences that may extend from $-\infty$ to $\infty$ reflected in its limits of the defining sum:
$$X(e^{j \omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} \tag{3}$$

When you want to make a calculation by hand, using symbolic expressions for signals, such as $x[n]=a^n u[n]$ and $h[n]=b^n u[n]$, you can calculate the results in time or frequency domains as outlined above.

Also, when you want to compute the same result by using a computer, you can use the time domain approach based on an LCCDE recursion (for IIR systems) or direct finite convolution sum (for FIR systems), BUT the frequency domain approach won't work with DTFT; as it's mainly a tool used to develop the mathematics of signals & systems theory, and it's not suitable for digital computer implementations, as its variable $\omega$ is a real continuous number.

What is used instead is the DFT (discrete Fourier transform) defined as $$ X[k] = X(e^{j\omega})|_{\omega = \frac {2\pi k}{N} } \tag{4}$$

where $k=0,1,...,N-1$ and $N$ is the length of the DFT, which's then called as an N-point DFT of signal $x[n]$.

Eq.4 implies that the DFT sequence $X[k]$ is obtained as the uniform samples of the DTFT $X(e^{j\omega})$, which is a periodic function, hence DFT $X[k]$ is also periodic but we only consider its first period from $k=0$ to $N-1$.

Since DFT sequences are inherently periodic, then their convolutions will also be periodic (circular). Therefore, whereas a linear convolution between aperiodic signals $x[n]$ and $y[n]$ is implied by the I-DTFT expression $$y[n] =\mathcal{I-DTFT} \{ X(e^{j\omega})H(e^{j\omega}) \}$$ instead a circular convolution between two periodic sequences is implied by the I-DFT expression $$\tilde{y}[n] = \mathcal{I-DFT} \{ X[k]H[k] \}$$

So, given that we want to compute a linear convolution between two aperiodic sequences $x[n]$ and $h[n]$ of lengths $L_x$ and $L_h$ respectively, using frequency domain by their $N$ point DFTs, $X[k]$ and $H[k]$, we actually have to compute a circular convolution between the periodic extensions of the signals $\tilde{x}[n]$ and $\tilde{h}[n]$ of periods $N$.

The key is in chosing a proper length $N$ of the DFTs, which must be long enough to avoid any time domain aliasing of the sequence $\tilde{y}[n]$, returned by the IDFT computation: $$ \tilde{y}[n] = \sum_{r=-\infty}^{\infty} y[n-rN] \tag{5}$$

where $y[n]$ is the result of the linear convolution that would be returned by the theoretical inverse DTFT and $\tilde{y}[n]$ is the periodic result of the periodic (circular) convolution implied by the inverse DFT.

Note that if any one of the signals are of infinite length, then it's NOT possible to compute their linear convolution using the DFT approach, as $N$ would go to infinity, practically impossible.
The implementation of a linear convolution via DFT then has the following steps:

Choose N according to the following criteria: $$N \geq L_x + L_h -1$$ which guarantees an alias-free reconstruction of the inverse signal $y[n]$ from its DFT $Y[k]$ of the computed circular convolution via $X[k]H[k]$.
Compute N-point DFTs $X[k]$ and $H[k]$ of $x[n]$ and $h[n]$.
Compute $Y[k] = X[k]H[k]$
Compute N-point inverse DFT of $Y[k]$ to produce the output $y[n]$


When you dont have the exact voice as a reference you can only go to estimate/measure the noise. Assuming you are facing white noise, the simplest thing is to look for either time-domain or frequency-domain parts of the signal, where no voice is present. In these parts there is only noise. Hence, if you measure the average signal power in these parts, you get the noise power P(noise). Then, you measure the overall signal power, which is P(voice) + P(noise). Finally, you can have
$$ SNR=\frac{P(voice)+P(noise)-P(noise)}{P(noise)} $$


Equalization and then Sharpening vs Sharpening and then Equalization

Since the Histogram Equalization isn't Linear Operator the order does matter and the result of the doing one before the other changes the final result.
Sharpening would work best when all the Dynamic Range is used.
Hence I'd first maximize the Dynamic Range by Histogram Equalization and then would apply the sharpening.
You might also think of the noise when making this choice.


Software implemented PLLs for signal synchronization, why use them?

For very high SNR signals locating the zero crossings will work pretty well.  You use PLL's when the signal doesn't have a high SNR, like the following.

As you can see by inspection, a zero-crossing algorithm wouldn't have a prayer of working on this signal.  A PLL, on the other hand, could do just fine.  That is, by the way, a tone with an SNR of 8 dB.


How do you calculate spectral flatness from an FFT?

If the definition of the flatness dictates that you use a power spectrum, then yes, you should square the magnitudes as the reference from the SciPy documentation indicates. In the equation that you referenced where you didn't see a squaring, I don't think you can read much into it; it says that

$$
S_{flatness} = \frac{\exp\left(\frac{1}{N} \sum_k \log (a_k)\right)}{\frac{1}{N} \sum_k a_k}
$$

but I don't see a definition for $a_k$ anywhere. If you want the spectrum to be proportional to the power in each bin, you need to square.


Is it possible to design a band-pass filter with fewer taps than a series HPF + LPF?

Assuming that your legacy FIR linear phase filters were designed using windowing method, then yes, you can combine those two FIR linear phase filters into a single FIR linear phase filter and it will be shorter than the the serial cascade.

However, you won't be convolving the individual impulse responses to get it. Rather, you should use a multi-band windowed FIR filter design technique. An example of how to do it can be found on Alan Oppenheim's DTSP book on chapter 7 filter design section.

Specifically he makes use of Kaiser window design approach acording to which the filter (window) length is empricially dependent on the transition band width and peak approximation error.

So, given the same transition band width and error ripples of a HPF -> LPF cacade, whose lengths are M samples each, you would be using a single multiband filter of length M with the same transition band and error ripple.

Note that as MattL also indicates, you may also chose different filter types to achieve similar respones, but that depends on what you want to achieve finally.


You are mixing two things modulation and multiple access. SC-FDMA has been chosen because it is battery-friendly modulation. SC-FDMA is used in both FDD and TDD modes.

There are several reasons why both FDD and TDD are defined.

Historically, people want to utilize efficiently the precious spectrum. FDD needs two bands while TDD needs only one.

Economically, LTE is a broadband technology, i.e. it is designed for high speed human data connection, and high speed human data connection favors downlink. For example in China, the dominant (and at the beginning of commercial LTE, the only) market of TDD LTE, Youtube video streaming takes a large part. In this case, reserving a band for UL is kind of wasteful. Moreover, LTE was commercialized in the middle of the Great Recession and operators were very skeptical in CAPEX spending.

Protocally, FDD and TDD are very different in term of multiplexing devices (UE), especially in HARQ-ACK reporting, UE scheduling etc. If you look at the procedure standards of LTE, just PHY layer for example, each mode requires very complicated scheduling. In mixing two modes in one cell, the level of complexity is huge for eNB, with no relevant benefit. However, mixing two modes between two or more cells are considerable and has been standardized since Release 11 in the Carrier Aggregation scope.

Technically, TDD supports better channel information (CSI) reporting thanks to the reciprocity of channel. Thus it is favored for beamforming, which requires accurate channel estimation, and a very good candidate for massive MIMO.

Finally, yes you can mix two modes if you want, but the efforts would not worth the benefit you get.


What about an algorithm in image domain?

For every pixel, check into the direction you desire to know (up-down in your images above, I guess) if there is a pixel with "the same" value (i.e. a value inside an acceptance band, of course. Let's say +/- 10%) as the pixel under investigation. If so, go further into that direction. Count every pixel that you find. If the number of counts exceeds a certain threshold (i.e. 1/4 FOV or so in your case), then there is a line.

I would first try to use something simple like that instead of using more complicated algorithms that might be suitable for different kinds of images. Yours pretty much look like MRI phase images that can be nasty ;) )


This is a rather general answer : I would say that, although each window function has specific characteristics that make it more or less suitable for specific situations, in this context you're not going to see a huge difference in the STFT of a windowed speech signal.

In particular, the two windows that you're comparing are quite similar, so I'd expect that you'd get pretty similar STFT results. You might try two windows that are quite different, like a rectangular and a Hamming or something like that, to see just how much difference those make in your STFT, and then use that difference as a sort of upper bound on the difference that you're likely to see from other windows.


SIFT: why it is invariant to affine tranlsation?

I consider that you already know the proof of how we get DoG from LoG, so I won't go further with it.

According to your question "how by approximation,using DoG, we get rid of that scale factor in SIFT". DoG works as a feature detector to detect blobs in the given image (blob like shapes) which you can say it is a bit kind of edge detectors here is a photo of some blobs detected using DoG:

These blobs detected at different scales using the scale space pyramid. Using larger values of sigma you get bigger blobs. Experimentally DoG give you the best notion of scale as you can see here :

So, experimentally DoG gives you the most robust features that are scale invariant by constructing the scale space and computing the DoG and finally extracting local maxima points.

Useful links:

1- Lecture PPT

2- Best known video lecture for SIFT

NOTE: This not a professional answer and I hope to hear more from other fellows


Do antisymmetric FIR filters preserve the shape of their input waveforms?

Such phase response is called Generalized linear phase, where you are allowed the $\pm \pi/2$ constant phase, and zero crossings (which add $\pi$ jumps in phase at one or more frequencies).

Filters with an antisymmetric impulse response do in fact have the $\pm \pi/2$ constant phase.

Simplest case: $h[n] = \delta[n] - \delta[n-1]$, with frequency response $H(e^{j\theta}) = 2j e^{-j \theta /2} \sin( \theta/2 )$.

In practical terms, these filters are types of differentiators: they all have a zero at $z=1$ (0 frequency) due to the fact that $\displaystyle\sum h[n] = 0$.

It is clear that these filters do not retain the waveform of the signal.

Normally you say that linear phase filters "retain the waveform" in case of low-pass filters, where you want to filter out noise but retain the shape of a pulse. More importantly, you want a good eye-pattern when decoding, and a non-linear-phase filter will result in asymmetric step responses and a poor eye pattern.


Periodicity of transfer function of FIR filter proof (Parks and Burrus, Digital Filter Design)


Understanding legacy code: Algorithm to remove radial lens distortion

Welcome to SP. I happen to have already answered your question in SO: see point 2 of my answer here.


IIR filter relationship between number of poles, order and second order sections

The filter order equals the number of poles, that's a fact. How people call and use parameters in software is a totally different question. Many routines that design standard IIR filters use frequency transformations to design bandpass or bandstop filters from a prototype lowpass filter. Often, the specified filter order is the order of the prototype lowpass filter, but the transformation doubles that order such that the resulting bandpass or bandstop filter has twice the order that was specified in the call of the routine. It works that way in Matlab/Octave. Maybe that's also the reason why you get twice the specified filter order with the SciPy routine and with the Iowa Hills software. Note that in both cases you obtain twice the specified filter order (or number of poles).


imfill would work with images and not volumes. However Matlab has many basic tools to fill the holes, that is to interpolate missing data points. What you are after is interpolation and Matlab has a function called interp3.


Question of Digital filter; stability and order

Some of your questions overlap so I’m not going to answer in the framework of 1 followed by 2, .....

An IIR filter ‘s stability is related to the location of its poles.  If the region of convergence contains the unit circle, the filter is stable. This definition covers both causal and anti causal filters, so you can have a stable causal filter and an unstable anti causal filter or vice versa.

Filters are implemented in terms of their coefficients so you need to root the denominator polynomial to find the poles and this is an ill conditioned problem.

see https://math.stackexchange.com/questions/70015/roots-of-bivariate-polynomials

for an example.  small perturbations in filter coefficients, like finite precision, can move the poles so that the region of convergence no longer contains the unit circle.

A good DSP book like Oppenhiem and Schaefer will discus work arounds like factorizing a polynomial into smaller stable polynomials and filter structures that are more robust to precision errors.

There are a number of functions in MATLAB that estimate oder for a given filter performance but these aren’t always reliable.

Ultimately one uses trial and error when designing a filter.  start with an oder and test the filter. if the filter meets spec, stop. if not, increase order by one and repeat.

How do you know if a filter is stable? you test it, which is what MATLAB is good at, although add ons like the fixed precision tool box might make things easier.

Overdetermined filters can be stable or unstable. An overdetermined filter can be a poor use of computation.

FIR filters are always stable. While an IIR filter requires fewer arithmetic operations than a FIR filter, some processors have execution pipelines that penalize recursion, so you can be doing more arithmetic faster with a FIR filter than the equivalent IIR filter.

A FIR filter can be linear phase, which can be desirable in many applications. An IIR filter can only be linear phase over a portion of bandwidth.

IIR filters can borrow from analog filters lke the Butterworth filter you mentioned

The answer to question 3 depends on what you mean by best.  It isn’t really an issue of order but how oder is structured


Matlab : How to plot the constellation diagram for this method?

The function you are using is from the Communication Systems Design Toolbox. You might have it not installed/licensed on your PC.

In case you just want to plot a constellation diagram, you can go with this code (Note, that it works for QAM constellations, not BPSK (since BPSK has only one domain)).

values = -7:2:7;  % Values in real and complex domain
[X, Y] = meshgrid(values, values);
QAM = X(:) + 1j*Y(:);
plot(QAM, 'x')

Actually, the QAM constellation contains all combinations of possible real and complex values. Hence, you need the cartesian product between real and complex parts. Taking this solution, you can plot the constellation diagram.


DFT calculations returning different results

Make sure you use the same scaling and X and Y axis so you compare apples to apples. The upper graph is in power the lower seems to be amplitude
Start with a known signal with a known answer. A sine wave or square wave for example and compare the outcome with the expected result.


Sum of Sine and Cosine with Random Phase as LTI System

So let's just look at it from a linear systems perspective:
$$
\begin{array}{ccl}
r(t) &=& Y_1(t) + Y_2(t)\\
&=& Z_1(t) \cos(2\pi f_0 t + \theta) + Z_2(t) \sin(2\pi f_0 t + \theta)\\
&=& [ X_1(t) * H_1(t) ] \cos(2\pi f_0 t + \theta) + [ X_2(t) * H_1(t)] \sin(2\pi f_0 t + \theta)
\end{array}
$$
where $X_1$ and $X_2$ are the respective (unlabeled) inputs, $H_1$ and $H_2$ are the impulse responses of the filters, and $*$ is convolution.

You appear to be asking whether the system is LTI when $H_1 = H_2 = H$ and $X_1 = X_2 = X$.  In that case, we get
$$
\begin{array}{ccl}
r(t) &=& [ X(t) * H(t) ] [\cos(2\pi f_0 t + \theta) + \sin(2\pi f_0 t + \theta)]\\
&=& \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)
\end{array}
$$

So it's linear:
$$
\begin{array}{ccl}
r_a(t) &=& \sqrt{2} [ a(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_b(t) &=& \sqrt{2} [ b(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_{a+b}(t) &=& \sqrt{2} [ (a(t) + b(t)) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4) = r_a(t) + r_b(t)\\
\end{array}
$$
But it's not time-invariant:
$$
\begin{array}{ccl}
r(t) &=& \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_\delta(t) &=& \sqrt{2} [ X(t-\delta) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r(t-\delta) &=& \sqrt{2} [ X(t-\delta) * H(t) ] \cos(2\pi f_0 (t-\delta) + \theta - \pi/4)\\
\end{array}
$$
as $r_\delta(t) \not = r(t-\delta)$.

EDIT

So, what happens if $X_1 = X_2 = X$ and $X\ \tilde{\ }\ N(0,\sigma^2)$ and $\theta\ \tilde{\ }\ U[0,2\pi)$?

I can see why $Y_1(t) = [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta)$ might be non-Gaussian.

However, I cannot see why $Y_1 + Y_2 = \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)$ will be Gaussian?

Can you please update your question to explain?


Getting blobs from similar colored contours in grayscale image

Here is a simple steps in c++ to identify the Objects you want to detect! Hope this solves your problem!

void main()
{
    Mat mSource_Bgr,mSource_Gray,mSmoothed,mEdge,mThres;
    mSource_Bgr= imread(FileName_S.c_str(),IMREAD_COLOR);
    cvtColor(mSource_Bgr,mSource_Gray,COLOR_BGR2GRAY);

    medianBlur(mSource_Gray,mSmoothed,5); /*Just to smooth the Edges*/
/*This is for Automatically calculating the Canny Threshold values*/
    double CannyAccThresh = threshold(mSmoothed,mThres,0,255,CV_THRESH_BINARY|CV_THRESH_OTSU);
    double CannyThresh = 0.1 * CannyAccThresh;

    Canny(mSmoothed,mEdge,CannyThresh,CannyAccThresh);

    Mat kernal_E,kernal_D;
    int erosion_size=2,dilation_size=2;
    kernal_E= getStructuringElement( MORPH_RECT,Size( 2*erosion_size + 1, 2*erosion_size+1 ),
                                       Point( erosion_size, erosion_size ) );
    kernal_D= getStructuringElement( MORPH_RECT,Size( 2*dilation_size + 1, 2*dilation_size+1 ),
                                       Point( dilation_size, dilation_size ) );
    dilate(mEdge,mEdge,kernal_D);//Just to make sure all the Edges are closed
    erode(mEdge,mEdge,kernal_E);

    imshow("mCanny_Edge",mEdge);
    imwrite(FileName_S+ "mCanny_Edge.bmp",mEdge);

    vector<vector<Point> > contours;
    vector<Vec4i> hierarchy;
    /// Find contours
    findContours( mEdge.clone(), contours, hierarchy, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_SIMPLE, Point(0, 0) );
    /// Approximate contours to polygons + get bounding rects and circles
    vector<vector<Point> > contours_poly( contours.size() );
    for( size_t i = 0; i < contours.size(); i++ )
    {
        approxPolyDP( Mat(contours[i]), contours_poly[i], 1, true );
    }

    for (int i = 0; i < contours_poly.size(); i++)
    {
        if(contours_poly[i].size()>100)//Just a filter You can add more here
            drawContours(mSource_Bgr,contours_poly,i,Scalar(0,255,0),2);
    }
    imshow("mResult",mSource_Bgr);

    }


Aperture Jitter in Digital to Analog Converters

Yes absolutely there is also the equivalent of aperture jitter results in a DAC. This is due to the effects of jitter in the sampling clock itself, and variation in the electronics of the time duration for the translation from the different digital levels to analog levels in the output.

In either case the jitter is not necessarily just due to a S/H circuit but the combined effects of that interface as well as the clock itself, and power supply noise can also translate to aperture jitter given the resulant changes in threshold levels.


What is the significance of first zero upcrossing of a cross correlation?

Not sure what you are asking specifically. The cross correlation of two cosine waves is just another cosine. The triangular envelope in the picture is simply a consequence of the limited length of the signals.

The phase of the cross-correlation cosine is given by the phase difference of the two input cosines. This phase difference determines the position of the zero crossings. So in effect the position of that zero crossing is some metric of the phase difference of the input signals


The IEEE used to host something called the Signal Processing Information Base, which is gone but there is a mirror at:

http://spib.linse.ufsc.br/database.html

There is also some links from Mike Porter's Ocean Acoustics Library at:

http://oalib.hlsresearch.com/


Alternatives to Hough transform for detecting a grid-like structure

An alternative to the Hough Transform would be the Radon Transform (1, 2). A rough description of an algorithm to detect a grid-like structure could look like this:

1. Perform Radon Transform from 0 to 180 degrees.
2. Find the two highest maxima in the angle bins.
3. For the two angles with maximal amplitude find the local maxima within the bin.
4. You can use the constraint that the maxima should have even spacing to deal with outliers.

EDIT:

Here is a small matlab snippet to illustrate step 1-3:

im = imread('grid.png');
[R, xp] = radon(im, 0:180);
imagesc(0:180, xp, R)

plot(max(R)) % the two maxima are at 65 & 117 degrees

plot(R(:, 65))

plot(R(:,117))

To answer your question from the comments: It appears to me from the one example image that you provided, that this method is more robust against small defects in the detected grid. Skeletons will rarely yield straight lines which might be a disadvantage for the subsequent Hough Transform.


Question regarding sampling rate of sine wave

The assumption here is the sine wave being samples is a known and assumed accurate reference, and the sampling rate is the unknown.
An algorithm that utilizes every sample with equal weight would provide for the best estimate assuming a white noise process. The approach the OP uses, and any similar "edge-detection" methods where we are making a decision on a cycle boundary and then measuring the distance between two such decisions, are prone to significantly higher errors under lower SNR conditions due to the sensitivity of noise on just the two decisions used; as done the noise effecting the result is the rms of the noise for each of the two decisions so is additive as well.
Such algorithms for determining the frequency of a single tone are well documented (see Question regarding estimation of signal tone parameters (frequency, amplitude, and phase) using Macleods algorithm). The approach would be to use these techniques with a scaled frequency designating the tone frequency to be determined as a frequency ratio relative to the sampleing rate $(f/f_s)$ rather than absolute frequency. Thus once $f/f_s$ is determined with such estimators, and with $f$ known, $f_s$ can be accurately determined.
The FFT is one such algorithm that uses every sample to provide a best estimate (under white noise conditions), but for the case of a single tone, the above references would provide more optimized solutions. See https://www.dsprelated.com/showarticle/1284.php where @CedronDawg worked through the math to derive the precise frequency estimate from two adjacent bins in the FFT result (for all the likely cases when the actual frequency is not centered on a bin), so that could be combined with the Goertzel algorithm which efficiently solves for a smaller subset of all DFT bins for a reasonable solution.