realtime software 90 degree phase shift a sinusoidal signal

If you have it in a buffer and have enough computing power, you can track it to quite significant accuracy using two DFT bins and the formulas shown here:

A Two Bin Solution

For minimum latency I would recommend a frame length near 2 and 1/2 cycles, though for a clean signal 1 and 1/2 cycles should also give good results.

As you do a series of sliding windows, this will generate a sequence of parameters.  A moving exponential average of the parameters should give you a good current value (when they are slowly varying) to use.  More sophisticated techniques are available.

Once you have the current parameters, you can reconstruct a signal with any phase lag and/or amplitude adjustment you want.

Or if you are tracking the phase value frame to frame, the most accurate fit is at the center of the frame ($t_f$ is at $n=N/2$).

$$ x(t) = A( t ) \cos( \Phi( t ) ) $$

Where $ \Phi( t ) $ is a (predictive) near linear function and $A( t_f )$ comes directly from the two bins solution.


Using continuous verses discrete wavelet transform in digital applications

As Mohammad stated already the terms Continuous Wavelet Transforms (CWT) and Discrete Wavelet Transforms (DWT) are a little bit misleading. They relate approximately as (Continuous) Fourier Transform (the math. integral transform) to DFT (Discrete Fourier Transform).

In order to understand the details it is good to see the historical context. The wavelet transform was originally introduced in geophysics by Morlet, and was basically a Gabor transform with a Window that grows and shrinks together with the selected scale/frequency. Later Daubchies (a physicist-ett from Belgium) realized that by choosing special orthogonal wavelet bases the infinitely redundant CWT can be critically sampled on a dyadic grid. From the resulting DWT the corresponding full CWT can be obtained by convolving the DWT with the reproducing kernel of the respective wavelet. The reproducing kernel is the CWT of the wavelet itself.

Daubchies findings gave a big boost to the wavelet theory in the early 80ies. The next big result was that the DWT can be computed very efficiently (this is sometimes called FWT [fast WT] as well) by using techniques from the theory of filterbanks, namely quadrature mirror filters (QMF) together with downsampling filterbanks. By constructing special QMFs the corresponding DWT can be computed via filtering and downsampling, which is the state-of-the-art algorithm to compute DWTs today. You do not need the scaling function to compute the DWT, it is just an implementation detail that FWT process.

Concerning the application side the CWT is the more ideal candidate for signal or time series analysis due to its more fine grained resolution and is usually chosen in most tasks (e.g. singularity detection). The DWT is more of interest in the context of fast non redundant transforms. The DWT has a very good energy compactification and is thus a good candidate for lossy compressions and signal transmissions.

Hope that clarified things.


A type 3 FIR filter has odd symmetry and an odd number of taps. For this reason it has frequency response zeros at $\omega=0$ (DC) and $\omega=\pi$ (Nyquist), corresponding to transfer function zeros at $z=1$ and $z=-1$:

$$H(1)=\sum_{n=-M}^{M}h[n]=0\\
H(-1)=\sum_{n=-M}^{M}(-1)^nh[n]=0$$

where $N=2M+1$ is the filter length. So it can neither be used as a high pass nor as a low pass filter.

Apart from that, the phase shift of $\pi/2$ caused by the odd symmetry is usually undesirable for frequency-selective filters. You could use such a filter for implementing Hilbert transformers or differentiators.


Design of a digital A-weighting filter with arbitrary sample rate

It's a common misconception that the approximation of an analog filter by a digital filter must be bad close to Nyquist. This idea might come from the ubiquity of the bilinear transform, for which this is usually indeed the case. Of course, there are certain constraints on the frequency response of discrete-time filters at Nyquist, but they do not necessarily need to result in a bad approximation of an analog filter in that frequency range. The quality of the approximation close to Nyquist depends on several factors, among which are the properties of the frequency response of the analog filter, and the fact whether only the magnitude or also the phase of the analog filter need to be approximated.

I've designed a 6th order IIR filter approximating the frequency response of an analog A-weighting filter, as defined here:

$$H(s)=\frac{k\cdot s^4}{(s+129.4)^2(s+676.7)(s+4636)(s+76655)^2}\tag{1}$$

with $k=7.39705×10^9$.

I chose a sampling frequency of $48$ kHz. The design procedure is a heuristic iterative procedure I came up with some time ago. It's a least squares approximation based on the equation error method, and I might write up all the details some day.

Below is a plot of the design result. Note that both plots go up to Nyquist ($24$ kHz). The left figure shows that one can't see any difference between the logarithmic plots of the analog and the digital frequency responses. The right-hand figure shows the approximation error, defined as the absolute value of the difference of the magnitude responses. The error shows a typical least squares behavior.

You can check out the filter yourself. Here are the coefficients:

b =

   0.169994948147430
   0.280415310498794
  -1.120574766348363
   0.131562559965936
   0.974153561246036
  -0.282740857326553
  -0.152810756202003

a =

   1.00000000000000000
  -2.12979364760736134
   0.42996125885751674
   1.62132698199721426
  -0.96669962900852902
   0.00121015844426781
   0.04400300696788968


Offset after Butterworth bandpass filtering

Your data have a non-zero DC component which is filtered out by the bandpass filter. So the filter output has approximately zero DC, which causes the corresponding vertical shift of the signal. If you want to retain the DC information in the signal you must use a lowpass filter instead of a bandpass filter.


How to Determine if FIR has a Linear_Phase Response w/o Matlab

Real-valued and causal FIR filters have a (generalized) linear phase response if and only if their coefficients satisfy either

$$h[n]=h[N-1-n]\tag{1}$$

or

$$h[n]=-h[N-1-n]\tag{2}$$

where $N$ is the filter length (number of taps).

There are four types of linear phase FIR filters, depending on the type of symmetry (even as in $(1)$ or odd as in $(2)$), and the filter length $N$ (even or odd). They have different properties, which is explained in this more detailed answer.


OpenCV/C++ connect nearby contours based on distance between them

If you are not worried about the speed or exact contour of hand, below is a simple solution.

The method is like this : You take each contour and find distance to other contours. If distance is less than 50, they are nearby and you put them together. If not, they are put as different.

So checking distance to each contour is a time consuming process. Takes a few seconds. So no way you can do it real time.

Also, to join contours, I put them in a single set and drew a convex hull for that set. So the result you are getting is actually a convex hull of hand, not real hand.

Below is my piece of code in OpenCV-Python. I haven't gone for any optimization, just wanted it to work, that's all. If it solves your problem, go for optimization.

import cv2
import numpy as np

def find_if_close(cnt1,cnt2):
    row1,row2 = cnt1.shape[0],cnt2.shape[0]
    for i in xrange(row1):
        for j in xrange(row2):
            dist = np.linalg.norm(cnt1[i]-cnt2[j])
            if abs(dist) < 50 :
                return True
            elif i==row1-1 and j==row2-1:
                return False

img = cv2.imread('dspcnt.jpg')
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
ret,thresh = cv2.threshold(gray,127,255,0)
contours,hier = cv2.findContours(thresh,cv2.RETR_EXTERNAL,2)

LENGTH = len(contours)
status = np.zeros((LENGTH,1))

for i,cnt1 in enumerate(contours):
    x = i
    if i != LENGTH-1:
        for j,cnt2 in enumerate(contours[i+1:]):
            x = x+1
            dist = find_if_close(cnt1,cnt2)
            if dist == True:
                val = min(status[i],status[x])
                status[x] = status[i] = val
            else:
                if status[x]==status[i]:
                    status[x] = i+1

unified = []
maximum = int(status.max())+1
for i in xrange(maximum):
    pos = np.where(status==i)[0]
    if pos.size != 0:
        cont = np.vstack(contours[i] for i in pos)
        hull = cv2.convexHull(cont)
        unified.append(hull)

cv2.drawContours(img,unified,-1,(0,255,0),2)
cv2.drawContours(thresh,unified,-1,255,-1)

Below are the results i got :


How to filter pulse like jumps in a signal?

I think that you will not be able to get acceptable results with linear filtering because the impulse noise covers the whole frequency range. You may want to look into nonlinear filtering techniques. As a first (and simple) approach try median filtering: http://en.wikipedia.org/wiki/Median_filter

I think if you optimize the window size of the median filter you should get better results than with purely linear filtering.


Block Processing (sometimes called Vector Processing) is a programming technique similar to, but not exactly the same as First-In-First-Out (FIFO) buffering.  It is double buffering of the samples coming into the process (be they from an A/D converter or a stream) and going out (to a D/A converter or another stream).  Block processing divides the processing time into slots of the entire block.  I.e. in audio, 32-sample blocks might be common in 48 kHz samplerate data.  So in one block time (which is 32/48000 second), the DSP algorithm is presented with 32 valid input samples at the beginning of the time period, and has the entire block time (of 32 sample periods) to compute 32 valid output samples.  This allows the cost of some of the work, like loading/saving states and loading coefficients and initializing accumulators and such, to be amortized over all 32 samples.  This makes the processing more efficient.  But the cost is 2 blocks of delay, due to the double-buffering.  And if there is feedback from one processing algorithm to another, that block of samples being fed back is from the previous block of samples, not the current block.  So feedback has a delay of 1 block of samples.
—
HRTF, I'll have to come back to that.  Maybe someone else can pick it up.
—
Sample rate is a well-known definition.  It is the number of samples per unit time expressed as frequency.  Usually it's samples-per-second expressed as Hz.  Bit depth is the same as the word-width of the samples.  CDs have 16-bit depth.  Decent audio processing DAWs have 24 or 32-bit depth.  Bit rate for uncompressed audio is simply the sample rate times the bit depth (times the number of channels).  It's the data rate of the audio stream.  For compressed audio, the bit rate is much reduces and continues to be the data rate of the audio stream.


First of all, I'm pretty sure that $H_e(k)$ is the $N/2$ point DHT of the even elements of the original sequence to be transformed, not of $H(k)$. The same is true for $H_o(k)$. As for the index $k$, when used with $H_e(k)$ and $H_o(k)$, it is evaluated modulo $N/2$, so it ranges from $0$ to $N/2-1$. So the index $k=N/2$ is mapped to the index $0$ when evaluating $H_e(k)$ and $H_o(k)$.


A large share of such questions are addressed via known sum of series. With trigonometrics arguments like $\cos(\alpha n+\beta)$, one solution is to use Euler/De Moivre formulae:

$$\cos(\alpha n+\beta) = \frac{e^{j(\alpha n+\beta)}+e^{-j(\alpha n+\beta)}}{2}$$

and recognize, under the Fourier sum, two geometric series:

$$\sum_{n=0}^{N-1} \frac{e^{j(\alpha n+\beta)}+e^{-j(\alpha n+\beta)}}{2}e^{-j2\pi kn/N}$$
with terms like :

$$\sum_{n=0}^{N-1}\left( e^{j(\pm\alpha -j2\pi k/N)}\right)^n$$
whose sum is well-known. For a given  $k$, you have:

$$
\begin{align}
\sum_{n=0}^{N-1} e^{\pm j(\alpha n+\beta)}e^{-j2\pi kn/N} &= & e^{\pm j\beta}\sum_{n=0}^{N-1} e^{\pm j\alpha n-j2\pi kn/N}\\
&= & e^{\pm j\beta}\sum_{n=0}^{N-1} e^{j(\pm \alpha -2\pi k/N)n}\\
&= & e^{\pm j\beta}\sum_{n=0}^{N-1} \left(e^{j(\pm \alpha -2\pi k/N)}\right)^n\\
\end{align}
$$
and then, noting

$$ \rho = e^{j(\pm \alpha -2\pi k/N)}$$
you can get:

$$\sum_{n=0}^{N-1}\rho^n =\begin{cases}N \textrm{ when } \rho=1, \\\frac{1-\rho^N}{1-\rho} \textrm{ when } \rho\neq 1.\end{cases} $$


Express Circular Convolution in Terms of Linear Convolution

$h[n]_N$ is the periodic extension of $h[n]$.

$$ h[n]_N = h[n + mN] \quad \text{for } 0 \le n+mN < N \text{ and } m \in \mathbb{Z}$$

so, given an $n$ and knowing $N$, you have to choose the correct integer $m$ so that

$$ 0 \le n+mN < N $$
or
$$ -n \le mN < N-n $$
or
$$ -\frac{n}{N} \le m < \frac{N-n}{N} = 1 - \frac{n}{N} $$.

Or, for the term $h[n-k]_N$,

$$ h[n-k]_N = h[n -k + mN] \quad \text{for } 0 \le n-k+mN < N \text{ and } m \in \mathbb{Z}$$

and

$$ -\frac{n-k}{N} \le m < 1 - \frac{n-k}{N} $$

so, you have to split your summation into two summations.


Frequency resolution of your spectrogram is connected to time window length and sampling frequency by following relation:

$$\Delta f = \dfrac{f_s}{N}$$

where $f_s$ is the sampling frequency, and $N$ is length of your time window in samples. $N$ can be easily expressed in seconds, in your case length of window is $T=0.5 \mathtt{s} $ and $f_s = 1000 \mathtt{Hz} $. Therefore window length in samples is: $$N=T\cdot f_s = 0.5 \cdot 1000=500$$
Which is giving you frequency resolution to be:

$$\Delta f = \dfrac{1000}{500}= 2 \mathtt{Hz}$$

You might also notice that substitution of $N$ with respect to $T$ and $f_s$ yields:

$$\Delta f = \dfrac{f_s}{T\cdot f_s}= \dfrac{1}{T} $$

where $T$ is obviously length of your window.


Comparison of Baseband simulation results and experimental results at Passband

Most of the time there is no good reason to simulate at passband, as long as you simulate the signal's and the channel's low-pass representation. Furthermore, assuming that you're using a modern radio such as an SDR, most of the hardware is operating at baseband anyway -- frequency up- and down-conversion are performed as close to the antenna as is feasible.

Also, be aware that, most likely, your simulation and experimental results will be hard to compare directly. For example, it is hard to find exactly what the SNR at the receiver is; you'll have to estimate it. Also, it's hard to determine the exact channel parameters. Other factors are non-linearities in the power amplifiers, synchronization imperfections, etcetera.


RNN as posterior probability estimation in speech recognition with HTK

It is very huge effort to build RNN-HMM system with HTK. It is easier to use more flexible toolkits like Kaldi, they have RNN-HMM already implemented.

HTK in latest version 3.5 beta already has support for DNN (not RNN but still a neural network), you can grep for HYBRIDHS in sources. For details on DNN-HMM it is better to read HTKBook 3.5 section 3.12 DNN-HMM Systems and Chapter 14 Artificial Neural Networks

You need to modify all places in HModel.c file to introduce your own type of the model. It is not a separate step you could introduce, you have to modify complex interaction between acoustic model and HMM search code.

It is also not straightforward to embed Python code into HTK C code, you might dump RNN scores in a file and read them from C, otherwise, you will have to reimplement your RNN code in C.


Is this statement correct from DSP aspect?

It's almost a matter of philosophy, i.e., difficult to argue hard facts.
On the one hand all the features you mention can be extracted from the raw signals. So in theory the network should be able to learn how to do that if they provide meaningful information for the task at hand. This is what part of the ML community is claiming: feature engineering is dead, let the networks figure out themselves what are the best features, we're only biasing them by what we think are good features, which may be wrong.
On the other hand, if you do have very good reasons to believe that some features are particularly useful, it can actually help the network. It wouldn't need to learn to discover such features first and thus, it could train much faster. This is related to the concept of known operator learning, where you embed fixed operators (that you know in advance should be beneficial to your task) into the network to make it train faster. I'm leaning towards this approach since it often makes sense.
If people used complex-valued STFTs to the input, no data were lost and all the features you mentioned (ZCR, F0, ...) could still be extracted from it. In a spectrogram, the phase is discarded, so one can argue that in fact some information is lost. If you just use your features as input though, it still does not represent the entire signal so some information is still lost and you may have to argue why you think this information is irrelevant.
Here is where it is difficult to argue hard facts and the only thing you may be able to do is show actual training results to make your point.


How to do a Fourier transform on input from a microphone and obtain the average amplitude of a certain frequency range

If you know the frequencies you're interested in, you might be better off just calculating the power / energy at those frequencies rather than using a full FFT.

The best algorithm to do that is probably the Goertzel algorithm.  There is a python implementation of that here, and a related SP.SE question here.


How Is Mixed Norm ($ {L}_{1, 2 }$) Better than $ {L}_{1} $ Norm for Sparse Representation?

The mixed norm allows you to impose some simple structure in the solution matrix. Using your example with $p=q=1$ then this means the solution could have arbitrary elements set to non-zero coefficients. This would not impose any structure in the solution. In this case the $L_{11} $ norm just sums the absolute value of all the matrix elements - similar to the Frobenius norm.

As an example, consider that $L_{11}$ norm of the following matrices are the same:

$\begin{bmatrix} 0&2 \\ 0 & 0\end{bmatrix}, \quad \begin{bmatrix} 0&0 \\ 2 & 0\end{bmatrix}, \quad \begin{bmatrix} 2&0 \\ 0 & 0\end{bmatrix}, \quad \begin{bmatrix} 0&0 \\ 0 & 2\end{bmatrix}$.

So there is no preference as to which element is set to non-zero

Instead, if $p=2$ and $q=1$, then the objective function is minimizing the $l_2$ norm of the columns and then the $l_1$ norm over the that. This means you're looking for an $a$ matrix with a sparse number of column vectors and where each column vector has a small $l_2$ norm. Of course this assumes you are using the $l_1$ norm as a proxy for sparsity.

Consider the $L_{2,1}$ and $L_{1,2}$ norms of following matrices:

$\begin{bmatrix} 2&0 \\ 2 & 0\end{bmatrix}, \quad L_{2,1} = \sqrt{8} \approx 2.8, \quad L_{1,2} = 4
$

and

$\begin{bmatrix} 2 & 0\\ 0 & 2 \end{bmatrix}, \quad L_{2,1}= 4, \quad L_{1,2} = \sqrt{8}
$

So minimizing the $L_{2,1}$ norm tends to prefer the structure of the first matrix over the second one, while minimizing the $L_{1,2}$ norm tends to prefer the structure second matrix over the first.

Thus the mixed norm allows you impose a soft constraint on the structure of the sparse solution that you're looking for. It's relatively easy to see that be you could also set it up to look for a sparse set of row vectors, each with a small $l_2$ norm. Deciding which to use is really dependent on how you set up your equations.


You are fitting a polynomial to a logarithmic decay curve, which is probably wrong.  You probably want to fit a log curve to a log curve.  You can use the naïve method of np.polyfit(log(x), y, 1) which will probably work fine, or you can weight the points to make it more accurate, or you can use scipy.optimize.curve_fit.  All are described here: https://stackoverflow.com/q/3433486/125507
The coefficients returned by the polyfit are the coefficients of the polynomial.  So with deg=2, you're fitting a parabola to the decay curve, which again is probably wrong: p(x) = p[0] * x**2 + p[1] * x + p[2]  With the np.polyfit(log(x), y, 1) method you'll be fitting a straight line to a straight line, so there will only be a slope and an offset: p(x) = p[0] * x + p[1]  The slope would then be proportional to the decay rate of the log curve (somehow, which you can figure out for the type of value you need).
After getting this working, you can probably vectorize it for speed, fitting all pixels simultaneously instead of using for loops.  np.polyfit can handle this, for instance: "Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column."


If you are using Python you could try this:

def S_point(signal, R_peaks):

    num_peak=R_peaks.shape[0]
    S_point_list=list()
    for index in range(num_peak):
        i=R_peaks[index]
        cnt=i
        if cnt+1>=signal.shape[0]:
            break
        while signal[cnt]>signal[cnt+1]:
            cnt+=1
            if cnt>=signal.shape[0]:
                break
        S_point_list.append(cnt)
    return np.asarray(S_point_list)

If you allready have the R-peaks then you can use the signal and the R-peaks to get the S_points with this function

S_points=S_point(signal,R_peaks)

And finaly plot it like this

plt.plot(signal)
plt.plot(R_peak,signal[R_peak], 'yo')
plt.plot(S_points,signal[S_points], 'ro')


Product of power signal and energy signal

If you are not looking for a rigorous proof, then the following would suffice for a particular class of power signals which has the property that their magnitude is bounded; i.e, $ |x_p(t)| \leq K$ for all $t$.

Numerous examples exist for such power signals like sinusoidal waves, complex exponentials, all kinds of periodic continuous-time signals which are defined via convergent Fourier series sums... Then for such a power signal $x_p(t)$ and a given energy signal $x_e(t)$ you can show that their product $y(t) = x_p(t) x_e(t)$ will be an energy signal as follows; given

$$ E_x = \int_{-\infty}^{\infty} |x_e(t)|^2 dt $$

and $ |x_p(t)| \leq K $ for all $t$ then the energy of y(t) is

$$ E_y = \int_{-\infty}^{\infty} |y(t)|^2 dt = \int_{-\infty}^{\infty} |x_p(t) x_e(t)|^2 dt \leq \int_{-\infty}^{\infty} K^2 |x_e(t)|^2 dt = K^2 E_x$$

Hence $y(t)$ will also be an energy signal.


Is this a correct interpretation of the DCT step in MFCC calculation?

Let me start from the beginning. The standard way of calculating cepstrum is following:

$$C(x(t))=\mathcal{F}^{-1}[\log(\mathcal{F}[x(t)])] $$

In the case of the MFCC coefficients case is a bit different, but still similar.

After pre-emphasis and windowing, you calculate the DFT of your signal and apply the filter bank of the overlapping triangular filters, separated in mel scale (although in some cases linear scale is better than mel):

In respect to cepstrum definition, you now represented the envelope of the spectrum (reduced spectrum) in mel-frequency scale. If you represent that, then you will see that it kind'a resembles your original signal spectrum.

Next step is to calculate the logarithm of the coefficients obtained above. This is due to the fact that cepstrum is supposed to be a homomorphic transformation that separates signal from the impulse response of the vocal tract, etc. How?

An original speech signal $s(t)$ is mostly convolved with an impulse response $h(t)$ of the vocal tract:

$$\hat s(t)=s(t)\star h(t)$$

In frequency domain convolution is a multiplication of spectra:

$$\hat S(f) = S(f)\cdot H(f) $$

That can be decomposed into two parts, based on following property: $\log(a\cdot b) = \log(a)+\log(b) $.

We also expect that impulse response is not changing over time, thus it can be easily removed by subtracting the mean. Now you see why we taking the logarithms of our band energies.

The last step in the cepstrum definition would be the Inverse Fourier Transform $ \mathcal{F}^{-1}$. The problem is that we have only our log-energies, no phase information, so after applying the ifft we get complex-valued coefficients - not very elegant for all this effort to be a compact representation. Although we can take the Discrete Cosine Transform, which is 'simplified' version of FT and get real-valued coefficients!
This procedure can be visualized as matching cosinusoids to our log-energy coefficients. You might remember that cepstrum is also called 'the spectrum of the spectrum'? That's the very step - we are searching for any periodicity in our log-energy envelope coefficients.

So now you see that now it's rather hard to understand how the original spectrum looked like. Additionally, we are usually taking only first 12 MFCC's, as higher ones are describing fast changes in log-energies, which is usually making recognition rate worse. So the reasons for doing DCT were the following:

Originally you must perform IFFT, but it's easier to get the real-valued coefficients from DCT. Additionally, we no longer have full spectrum (all frequency bins), but energy coefficients within mel filter-banks, therefore usage of IFFT is a bit of overkill.
You see on the first figure that filter banks are overlapping, so the energy from ones next to each other is being spread between two - DCT allows to decorrelate them. Remember that this is a good property for example in case of Gaussian Mixture Models, where you can use diagonal covariance matrices (no correlation between other coefficients), instead of full ones (all coefficients are correlated) - this simplifies things a lot.
Another way of decorrelating mel frequency coefficients would be PCA (Principal Component Analysis), technique solely used for this purpose. For our luck it was proven that DCT is a very good approximation of PCA when it comes to decorrelating signals, hence another advantage of using Discrete Cosine Transform.

Some literature:

  Hyoung-Gook Kim, Nicolas Moreau, Thomas Sikora - MPEG-7 Audio and
  Beyond: Audio Content Indexing and Retrieval


Correlation: A Magnitude or Power quantity?

There are different ways to look at this.
In general, it's useful to think of the correlation as a power quantity.

If you slog through the units, the correlation will have power-like units: $V^2$ for example
The Fourier Transform of the autocorrelation is the Power Spectral Density, which is clearly a power-like quantity

So what's happening in your example? Let's assume that we have $y = x + n$ where $x$ and $n$ are mean-free, uncorrelated and have the variances $P_x = \sigma_x^2 = <x^2>$ and $P_n = \sigma_n^2 = <n^2>$. Since they are uncorrelated we have $P_y = P_x+P_n $
The SNR would be $$SNR = 10log_{10}\frac{P_x }{P_n}$$
The argument to the log is a ratio of powers.
The covariance is
$cov(x,y) = <x\cdot (x+n)> = P_x$ which has power-like units. The normalized correlation coefficient comes to be.
$$\rho = \frac{cov(x,y)}{\sqrt{{P_x \cdot P_y }} = \frac{P_x}{\sqrt{P_x \cdot (P_x + P_n) }} $$
$$\rho = \frac{cov(x,y)}{\sqrt{P_x \cdot P_y }} = \frac{P_x}{\sqrt{P_x \cdot (P_x + P_n) }} $$
So that's also a ratio of powers: cross-power to geometric mean of the signal powers. Since it's a ratio it's unitless. Let's square it:
$$\rho^2 =\frac{P_x^2}{P_x \cdot (P_x + P_n) } = \frac{P_x}{P_x + P_n } $$
The most reasonable physical interpretation of this is STILL a ratio powers: signal power to signal power plus noise power. Since the ratio is unitless, the square of ratio is still unitless, so there is no direct conflict here.


Analytically designing a notch-filter for specified frequency 50 Hz

You've designed a causal filter with a notch at $\omega_0=100\pi$. But the result is probably not what you want. Note that you've designed an FIR (finite impulse response) filter. Its frequency response has a large overshoot towards high frequencies.

What you actually want is an IIR filter, with poles away from the origin of the complex $z$-plane. A simple way to design a notch filter is to place the zeros at the desired notch frequency (as you did), and place poles at the same angle but slightly inside the unit circle. If $0<r<1$ is the pole radius, your transfer function looks like this:

$$H(z)=\frac{(z-e^{j\omega_0})(z-e^{-j\omega_0})}{(z-re^{j\omega_0})(z-re^{-j\omega_0})}=\frac{z^2-2z\cos(\omega_0)+1}{z^2-2zr\cos(\omega_0)+r^2}\tag{1}$$

You can apply scaling to the transfer function $(1)$ such that its gain at DC equals 1.

Finally, note that you have to take the sampling frequency into account. The notch frequency $\omega_0$ in $(1)$ is normalized. If $f_0$ is the desired notch frequency in Hertz, and if $f_s$ is the sampling frequency in Hertz, then

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

As an example, take $f_0=50\text{ Hz}$, $f_s=1000\text{ Hz}$, and $r=0.98$. The figure below shows the transfer function of the FIR notch filter you proposed (top), and the IIR transfer function according to Equation $(1)$ (bottom):


What is the minimum sampling rate of a $\sin(\pi t)$ signal?

$$F_\text{s} = 1 + \epsilon$$

where $\epsilon > 0$.  you must sample at some frequency exceeding twice the bandwidth.  (all this assmuming real $x(t)$ and none of this bandpass sampling stuff.)

  When the sampling frequncy is twice that of the frequency, the spectrum is not present and in time domain all the samples we get is 0.

nooooooo! that's not the case!  you could miss-sample it by $90^\circ $ and

$$x[n] = x(nT) = \cos(\pi n) = (-1)^n$$

instead.

  In this case does the sampling rate should be greater than $2 f_\text{s}$? then how does Nyquist rate fit in this

two things:

the original Nyquist paper (i'll try to find a link to it), actually excludes pure sinusoids, particularly a pure sinusoid sitting at Nyquist.
to deal with the $x[n] = x(nT) = \cos(\pi n + \phi)$ thing (a sinusoid at the Nyquist frequency), you have to accept the notion of Dirac delta functions, $\delta(f)$ in the frequency domain.  and then you will see that, at Nyquist, ($f = \frac{F_\text{s}}{2}$), there will be overlapping and aliasing.  but unlike other aliasing, you lose only half of your information, just the phase.  so if there were reconstruction done, it would be only on the basis of the quadrature component of the

$$ x(t) = cos(\pi t + \phi) = \cos(\phi)\cos(\pi t) - \sin(\phi)\sin(\pi t) $$

input. if $F_\text{s} \triangleq \frac{1}{T} = 1$ and $x[n]=x(nT)$, only the $\cos(\phi)\cos(\pi t)$ term counts.


Looking for project ideas to match my interests

I am a former physician, having done half graduation in mechanical engineering, then a master degree in product design. Nowadays, I work for a biomedical/bioelectric-equipment company, doing mostly software design (algorithms and workflow), using Numpy (matlab like) for prototyping.

THAT fulfills my need (incredibly similar to yours) to solve new problems, work with pretty-picture generation via numbercrunching, programming, and putting the head to perform hard intellectual work very often. I can tell no two months from the three years I am working here have been like one another, and the product portfolio keeps increasing.

In my spare-time, I have interest in geospatial stuff, most probably because I like to ride bikes on relatively remote places with maps and GPS, and the very description of geographical stuff allows for a lot of pretty-picturing. Recently, the avalanche of bike-riding GPS data on sites such as Strava, and the current problems CAUSED by traffic-engineering in big cities makes me think these statistical analysis from people displacement's big data can form a body of knowledge for an "anti-traffic-engineering" oriented for actual people instead of private automobiles.

I can "feel" there's a huge future demand, arriving pretty soon, for software tools that transform datasets into useful information. That can happen with medical imaging, geo-statistical imaging, traffic-engineering micro-modelling, and a lot of fields where data collection, transformation, analysis and interpretation is of huge importance, but still needs further development.

I'm not sure I am actually answering your question, but I feel your question is open to some daydreaming, so I hope this helps, and would be glad to answer any comment you may have!


Motivation for using gabor transform instead of fourier transform

The fourier transform gives you very fine resolution in the frequency domain, but during the transformation, you loose all the information about when (for time signals) or where (for images) these frequencies occur in your input signal.

The Gabor transform alleviates this problem by windowing the base functions of the fourier transform with a Gaussian impulse. The transform gives the similarity between the signal and the time- and frequency-shifted window, so you may get a sense when/where certain frequencies occur in your input signal.

However, the windowing comes at the cost of losing some of the frequency resolution. The Gabor  transform is optimum in the sense that it's windowing function yields the minimum product between time and frequency resolution.


How can I move a signal in frequency so that it ends at the frequency zero?

One way of moving a frequency response of signal is multiplying it with a sinusoid. Assume x(t) is your signal.Let the sinusoid be $$cos(w_0t)$$
If I multiply this to x(t) as $$y(t)=x(t)cos(w_0t)$$
Then Fourier transform of y(t) is $$Y(w)=0.5(X(w-w_0)+X(w+w_0))$$
This is because multiplying in time domain of two signals is equivalent to its convolution in frequency domain.Since fourier transform of the above sinusoid is$$C(w)=\frac{\delta(w-w_0)+\delta(w+w_0)}{2}$$Due to the sifting property of delta function we get the above Y(w) after multiplying.
Although here we get two signal the required one can be obtained by passing it through band pass filter


The conventional definition of the DFT for a length $N$ signal (without zero-padding) is

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$

So there is no scaling involved. Scaling is applied to the inverse DFT:

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

If you pad a length $N$ signal with $M$ zeros, the resulting DFT length is increased to $N+M$, and the DFT becomes

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/(N+M)}\tag{3}$$

The corresponding inverse DFT is

$$x[n]=\frac{1}{N+M}\sum_{k=0}^{N+M-1}X[k]e^{j2\pi nk/(N+M)}\tag{4}$$

EDIT:

In reaction to your comments, it is indeed necessary to sum over all indices $k$ up to $N+M-1$ in the sum of Eq. $(4)$. A simple example will illustrate this. Take a sinusoidal signal with an integer number of periods inside the DFT length:

N = 16;
n = 0:N-1;
x = sin(pi/4*n);
X = fft(x);

Its DFT X consists of two impulses. However, the DFT of the zero padded signal shows the underlying sinc-like behavior due to the rectangular window. Note that the DFT is just a sampled version of the DTFT, which has a continuous frequency variable. Zero padding just increases the sample density in the frequency domain:

xz = [x,zeros(1,N)];
Xz = fft(xz);
subplot(1,2,1),stem(abs(X))
subplot(1,2,2),stem(abs(Xz))


Is $O (N \log N)$ FFT speed the fastest we can ever attain?

While the general problem remains an open problem in theoretical computer science, there are (randomized) algorithms known as sparse fast fourier transforms which can approximate DFTs to a certain approximation factor when there are only a few DFT coefficients you are concerned with.

You can find some details in these papers and related papers (downloadable here):

Hassanieh, Haitham, et al. "Nearly optimal sparse Fourier transform." Proceedings of the 44th symposium on Theory of Computing. ACM, 2012.
Hassanieh, Haitham, et al. "Simple and practical algorithm for sparse Fourier transform." Proceedings of the Twenty-Third Annual ACM-SIAM Symposium on Discrete Algorithms. SIAM, 2012.


Cost function for LTI system identification

The chosen cost function is the mean squared error, i.e., the integral over a squared magnitude of the difference between frequency responses. The function

$$E(e^{j\omega})=H(e^{j\omega})-\frac{B(e^{j\omega})}{A(e^{j \omega})}\tag{1}$$

depends on frequency, so you can't minimize it directly, unless you want to minimize it for exactly one frequency $\omega$, which is of course pointless. You can choose several error measures depending on $E(e^{j\omega})$ given in $(1)$. Two common choices are

$$\varepsilon_1=\max_{\omega}W(\omega)|E(e^{j\omega})|\tag{2}$$

and

$$\varepsilon_2=\int_{0}^{\pi}W(\omega)|E(e^{j\omega})|^2d\omega\tag{3}$$

with some positive weighting function $W(\omega)$. Note that $\varepsilon_1$ and $\varepsilon_2$ given by $(2)$ and $(3)$ do not depend on frequency.

The authors of the paper you refer to chose to minimize the weighted mean square error given by $(3)$. However, instead of using the linear difference $(1)$, they chose to minimize the average logarithmic difference, i.e., the average error on a dB scale.

The problem with the minimization of $(3)$ is that for IIR filters, it results in a non-linear optimization problem, which is much harder to solve directly, and which might also have locally optimal solutions that are far from the global optimum. The cost function $\hat{C}$ in your question is linear in the filter coefficients. The point of the iteration is now to solve a sequence of linear minimization problems (which is simple, just solve a system of linear equations) in order to compute the solution of the originally non-linear optimization problem (the minimization of $(3)$). Note that if convergence is achieved, then $A_{i-1}(e^{j\omega})=A_{i}(e^{j\omega})$, so the cost function $\hat{C}_i$ is identical to the cost function of the original non-linear problem. Yet, only linear minimization problems are solved in each iteration.

I'm not sure I completely understand your last question, but the weighting function is there to change to problem from minimizing the average squared difference between frequency responses to the average squared differences between the logarithms of the magnitude responses. The necessary weighting function is unknown, but if the procedure converges - note that there is no guarantee that it does in all cases - the final weight function is such that the mean squared logarithmic error is minimized. Note that this latter problem is highly non-linear, whereas the proposed procedure only solves linear subproblems.


Practical Applications of upsampling and downsampling

Simplest reason for downsampling is simply to reduce the amount of data you have. Say you have a audio raw stream sampled at 44.1 kHz. You may simply want to reduce your data rate to 22.05 kHz because your algorithm can work with that sampling rate. Apply the same principle anywhere else. Just because you have more data does not mean you always need it. It may sound untrue but it is like saying do I need something upto 10 decimal points of accuracy or will 2 decimal points do? [Not exactly the same but you get the idea]

A simple reason to upsample to to interpolate. You have data which is sampled at a particular rate. You want to upsample it to get intermediate points between the actual data.  [Think of curve fitting. you want more points on the curve.]

There are implications to downsampling and upsampling. Not every downsampling is the same and not every upsampling give accurate interpolation. You will learn more as you go along.

Edit: Thought I'd give another simpler real life example
Image scaling: if you have a small image. You upscale it. This is upsampling.
If you want a smaller version you resize to a smaller size. This is downsampling.

Ofcourse practical scalers will not just add 0's and upsample and downsample. [You could do that too, but the output will not look good]. But it should give you another idea of where upsampling and downsampling is done.


Bessel filter second-order sections Q and Fc multiplier derivation

There is no simple closed-form equation to calculate Bessel filter bi-quad coefficients. The poles of the filter come from a Bessel polynomial. Higher order Bessel polynomials are determined using a recursion relationship. You need computer algorithms to calculate the poles and resolve the bi-quad coefficients.

Two books I can recommend: “Passive and Active Filters” by Wai-Kai Chen and “Analog Filter Design” by M.E. Van Valkenburg.


Protocols like w-mbus or Zigbee typically use some form of frequency-spreading like DSSS. In DSSS we do not send symbols directly using some pulse shape and encoding the information in the sign but instead we use a rapidly alternating sequence of $\pm 1$ as a pulse. This leads to an increase in the occupied bandwidth ("spreading"). The receiver multiplies the incoming signal with the same sequence, coherently combining back all the information ("de-spreading"), which overall increases the robustness to fading and reduces the peak power spectral density. See the link for more details.

What it means practically is that each symbol is further sub-divided in time domain into smaller units called "chips". The rate at which these chips are sent is called chip rate which itself is a measure for the occupied bandwidth. For instance, a chip rate of 500 kchips/s would mean that 500'000 chips per second are sent, which requires something like 500 kHz (practically more, due to the required guard bands).

If you want to know your data rate, you need to know the code (i.e., how many chips per symbol are used) and the modulation and coding scheme (i.e., how many bits per symbol are sent).

For example, this TI reference design states 100kcps with 11-bit DSSS coding and 2-FSK modulation. This gives 100kchips/s / 11chips/symbol * 1bit/symbol ~ 9kbit/s (minus packet overhead).


Both are orthogonal in your example, in different ways. The FDM example provides orthogonality by properly spacing the signals in the frequency domain (their dot product is 0). OFDM is far more spectrally efficient where you have overlapping carriers that are orthogonal.

From an implementation standpoint OFDM can be achieved with a DFT. Each of the carriers is a much lower bit rate that turns a frequency selective fading channel into flat fading. This is not true for your FDM example. There are many other differences but these are a few.


Segmentation and Labeling of White Blood Cells

I can't recommend anything in the realm of your clustering algorithm, but this looks to be pretty simple. The initial segmentation worked nicely; a simple watershed should separate those cells.

For example, in MATLAB:

function W = dowatershed(I)
    %% I is a logical mask (i.e., im2bw)
    d = -bwdist(~I);
    l = watershed(d);
    W = I;
    W(l == 0) = 0;
end


There is both a hue shift and a saturation shift from the full color image to the band artifact. Converting the color image to grayscale generates an image that still has an artifact, but that is less obtrusive.

There are "image inpainting" algorithms that you might try to correct the artifact, though it would help to first isolate the affected region. Here are a few thoughts on identifying the artifact:

Run Hough or RANSAC line fits for horizontal and vertical lines in the image. For a natural scene, and especially in a scene of the Earth, it would be odd to have horizontal or vertical lines corresponding to real features. Even the structures visible from space such as canals and large constructions tend to meander.
Try a form of image subtraction that subtracts a grayscale image from the original color image. Since the artifact is more grayscale-like and desaturated, you should be able to detect a region of hue and saturation differences. I would recommend working in HSV space for this operation. Run an algorithm in the vertical direction to look for a peak or valley in the resulting (color HSV - gray HSV) image.
Test row-by-row cumulative changes in Hue and Saturation. From row N to row N+1, the same X pixel generally won't change much in H, S, or V. However, you should notice a significant difference between the row just before the artifact and the row at the artifact. The difference might be as simple as the root mean square of H and S differences at each pixel, or even just the sum of differences. Or whatever. Pick a metric that is simple and that works.
Double check the original, raw image to determine whether the horizontal artifact extends into the black region of space on the left side of the image. In the copy of the image I downloaded, it appeared that the artifact did not extend into the black region.

Here are specs for JunoCam:
http://en.wikipedia.org/wiki/JunoCam

"The camera uses a Kodak image sensor, the KODAK KAI-2020, capable of color imaging at 1600 x 1200 pixels. It has a field of view of 18 x 3.4 degrees with three filters to provide color imaging."

That sensor should give a nice color image without weird artifacts when used on earth in gentler circumstances.

Although the artifact did not appear to extend into the black region of space at left, that could be a matter of dynamic range/windowing: there may be some difference in the charge on those pixels, but perhaps that difference is clipped to zero.

My guess is that this is a readout problem as the image rows are dumped. With older analog cameras we could see all sorts of crazy artifacts related to EM interference. With digital cameras that's usually much less of a concern, but one possibility is that some short, localized EM event occurred that somehow affected just certain rows.

If you zoom in to the top of the band artifact, you'll see desaturated rainbow artifacts. There are two non-adjacent rows that are brighter than their neighbors, too.

Since the artifact is lined up so nicely with the horizontal rows of the image, I would expect the problem is some internal electronics problem rather than a problem with lensing, color filters, flaring, or anything like that. Unless the camera became damaged by a particularly strong EM event, I would guess that the artifact would occur only rarely.


Why is the Fourier (or cosine) transform decorrelating?

Whether that decorrelation happens depends on the signals you put in – it's not a general property of the DFT!

Especially, when you model your signal as sum of narrowband signals, then you'll find that due to the DFT of a narrowband signal being concentrated on one or very few bins, and practically 0 elsewhere, the decorrelation simply stems from the fact that it's easy to represent the same signal in frequency domain with fewer non-zero coefficients than in time domain.

The fact that these non-zero coefficients are then uncorrelated amongst each other depends again on the signal model, and isn't widely true – for example, if you have a line spectrum DFT, for example from observing a rectangular wave digital signal, the lines aren't uncorrelated at all.


Real Time Object Tracking using image processing

What kind of an object?  If it is a compact object moving in some predictable fashion, you can use the Kalman filter to track based on its motion.  If the object has some texture, and doesn't move very fast relative to the frame rate, you can use the KLT (Kanade, Lucas, Tomasi) algorithm to track points on it.  If it has a distinctive color pattern you can use a Mean-shift tracker to track it.


Differences in transforming a sinc FIR filter from lowpass to highpass using spectral inversion and spectral reversal

This answer is gonna come in installments.
First consider the ideal $\operatorname{sinc}(\cdot)$ filter, which in the frequency domain is an ideal brickwall.  We'll do this in normalized angular frequency, $\omega$.
Let
$$ H_\mathrm{LPF}(e^{j\omega}) = \begin{cases}
1, \qquad  |\omega| < \omega_0 \\
0, \qquad  \omega_0 <  |\omega| \le \pi \\
\end{cases} $$
and define the frequency response outside of $\pm\pi$ to be periodic:
$$ H_\mathrm{LPF}(e^{j(\omega+2\pi)}) = H_\mathrm{LPF}(e^{j\omega}) \qquad \forall \omega \in \mathbb{R} $$
This will have impulse response of
$$ h_\mathrm{LPF}[n] = \frac{\omega_0}{\pi} \operatorname{sinc}\left( \frac{\omega_0}{\pi} n \right) $$
Now, if you make a high-pass filter by subtracting the low-pass from a wire:
$$ H_\mathrm{HPF}(e^{j\omega}) = 1 - H_\mathrm{LPF}(e^{j\omega}) $$
Then the cutoff frequency remains the same $\omega_0$ and the impulse response is
$$\begin{align}
 h_\mathrm{HPF}[n] &= \delta[n] - h_\mathrm{LPF}[n] \\
 & = \delta[n] - \frac{\omega_0}{2} \operatorname{sinc}\left( \frac{\omega_0}{\pi} n \right) \\
\end{align} $$
Now, suppose you make a high-pass filter by heterodyning the low-pass from DC to Nyquist:
$$ H_\mathrm{HPF}(e^{j\omega}) = H_\mathrm{LPF}(e^{j(\pi+\omega)}) $$
Then the cutoff frequency becomes $\pi-\omega_0$ and the impulse response is
$$\begin{align}
 h_\mathrm{HPF}[n] &= (-1)^n h_\mathrm{LPF}[n] \\
 & = (-1)^n \frac{\omega_0}{\pi} \operatorname{sinc}\left( \frac{\omega_0}{\pi} n \right) \\
\end{align} $$
So, if we harmonize the two definitions of the cutoff frequencies, the two brick-wall high-pass filters will be identical in the frequency domain, so they should in the time domain.  But that would mean that
$$\begin{align}
\delta[n] - \frac{\omega_0}{\pi} \operatorname{sinc}\left( \frac{\omega_0}{\pi} n \right) &\stackrel{?}{=} (-1)^n \frac{\pi - \omega_0}{\pi} \operatorname{sinc}\left( \frac{\pi-\omega_0}{\pi} n \right) \\
\\
&= (-1)^n \frac{\pi - \omega_0}{\pi} \frac{\sin\left( \pi \frac{\pi-\omega_0}{\pi} n \right)}{\pi \frac{\pi-\omega_0}{\pi} n} \qquad & (n \ne 0) \\
\\
&= (-1)^n \ \frac{\sin\big( (\pi - \omega_0) n \big)}{\pi n} \qquad & (n \ne 0) \\
\\
&= (-1)^n \ \frac{\sin(\pi n)\cos(\omega_0 n) - \cos(\pi n)\sin(\omega_0 n)}{\pi n} \qquad & (n \ne 0) \\
\\
&= (-1)^n \ \frac{0\cos(\omega_0 n) - (-1)^n\sin(\omega_0 n)}{\pi n} \qquad & (n \ne 0) \\
\\
&= (-1)^n \ \frac{- (-1)^n \ \sin\left(\pi\frac{\omega_0}{\pi} n\right)}{\pi n} \qquad & (n \ne 0) \\
\\
&= \frac{- \sin\left(\pi\frac{\omega_0}{\pi} n\right)}{\pi n} \qquad & (n \ne 0) \\
\\
&= -\frac{\omega_0}{\pi} \frac{\sin\left(\pi\frac{\omega_0}{\pi} n\right)}{\pi \frac{\omega_0}{\pi} n} \qquad & (n \ne 0) \\
\\
&= -\frac{\omega_0}{\pi} \operatorname{sinc}\left(\frac{\omega_0}{\pi} n \right) \qquad & (n \ne 0) \\
\end{align}$$
So I guess it's true, because you can also show that equally exists for the case of $n=0$.


If this represents the received signal where $a$ and $b$ are the real and imaginary parts respectively $s(t-\Delta) \approx [a \cos{(2\pi f_o \Delta)}-b \sin{(2 \pi f_o \Delta)}] \phi(t) \sqrt{2} \cos{(2 \pi f_o t)} + [b \cos{(2 \pi f_o \Delta)} + a \sin{(2 \pi f_o \Delta)}] \phi(t) \sqrt{2} \sin{(2 \pi f_o t)}$

one can see exactly how the signal has been rotated. Therefore the rotation can be countered by introducing a sine or cosine to counteract the delay.


Taking audio frequency and picking out peaks (Peak Detection)

I found the answer finally. I found a great article that explains many different libraries that can be utilized for peak detection.  I now have the peaks I am really interested in, and can now create the output I am requiring.

Finding Peaks in Python


Complex conjugation in the discrete Fourier domain should result in spatial domain horizontal and vertical mirroring around the origin. The origin happens to be at row 1 column 1, as is evident from MATLAB's definition of FFT, with Y = fft(X):

$$Y(k)= \sum^n_{j=1} X(j) e^{−2πi(j−1)(k−1)/n}$$

For $j = 1$ the phase is zero for all $k$ so anything on column 1 or on row 1 will stay there through complex conjugation. Column 11 and row 11 are each 10 steps from the origin horizontally or vertically, as are column 12 and row 12 if you measure the distance across the border. So the reversal swaps the data between columns 11 and 12, and between rows 11 and 12. Figure 1 below illustrates the mirroring.

Figure 1. Row or column indices and the corresponding phases of frequency $k = 2$ plotted as angles on the circle.

If the dimensions of your data were even rather than odd as they are now, there would be a middle row and a middle column that would be at equal distance from the origin both directly and across the border.


How to compare 2 MFCC (compare them if they where created from same speaker)

You can find the distance between the MFCC features for different speakers. The distance should be zero for the same speaker and larger with different speakers.
You can use different types of distance functions or even similarity functions used generally to compare between probability density functions.
Also, corresponding to each speaker you should generate a single matrix of MFCC coefficients which is composed of multiple coefficients for each frame. So, the output MFCC should be of dimension: (number of frames, number of coefficients).
Every speaker's MFCC dimension should be the same to compare them.


The Gradient of Least Squares of 2D Image Convolution

The easiest approach would be writing each case using Matrix Form of the convolution.
In this answer we assume the discrete convolution is applied only on valid support (Matching MATLAB's valid parameter for the convolution).
Namely, given $ x \in \mathbb{R}^{m \times n} $ and $ h \in \mathbb{R}^{k \times l} $ then $ h \ast x \in \mathbb{R}^{ \left( m - k + 1 \right) \times \left( n - l + 1 \right) } $. Needless to say $ m \geq k $ and $ n \geq l $ as otherwise the operation isn't well defined. Pay attention that this form of convolution isn't commutative.

Remark
For full convolution the problem is easier as the operation is commutative hence no difference between the gradient with respect to $ h $ or $ x $ as one could switch them.

Gradient with Respect to Convolution Kernel $ h $

The matrix form is given by:

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

Where $ X $ is the 2D Convolution Matrix Form of the image. Then:

$$ \frac{\mathrm{d} f \left( h \right) }{\mathrm{d} h} = {X}^{T} \left( X h - y \right) $$

The $ {X}^{T} $ forms a correlation (Versus Convolution) with full support of the operation (Equivalent of the full convolution shape in MATLAB syntax).

Hence we have:

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

In MATLAB Code:

conv2(mX(end:-1:1, end:-1:1), (conv2(mX, mH, CONVOLUTION_MODE_VALID) - mY), CONVOLUTION_MODE_VALID);

Gradient with Respect to Convolution Image $ x $

The matrix form is given by:

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

Where $ H $ is the 2D Convolution Matrix Form of the kernel. Then:

$$ \frac{\mathrm{d} f \left( x \right) }{\mathrm{d} x} = {H}^{T} \left( H x - y \right) $$

In MATLAB Code:

conv2((conv2(mX, mH, CONVOLUTION_MODE_VALID) - mY), mH(end:-1:1, end:-1:1), CONVOLUTION_MODE_FULL);

MATLAB Code

The full code is available on my StackExchange Signal Processing Q59089 GitHub Repository.


FFT calculation suitable, Wav File, Java

To show only one spectrum then FFT of the whole sample could be done.But for better result you could take a window and take FFT. For taking windows of 4096 with 50% overlap then it is equivalent to doing FFT
no of times
$$
\frac{480000}{4096*0.5}
$$
times that is equivalent to
$$
\frac{Total\hspace{1mm} sample\hspace{1mm} length}{(window \hspace{1mm}length)*(amount\hspace{1mm} of\hspace{1mm} overlap)}
$$


Intuitively, what is fourier series representation of a signal? Also intuitively what is frequency response?

A fourier series essentially breaks apart a periodic signal to represent it as an infinite sum of sine waves that are in that signal. So a fourier transform of a signal breaks apart the signal to look at all the frequencies that are in that signal (essentially looking at all the ingredients that make up that signal, ingredients being the frequencies).

A frequency response is the frequencies of an output signal compared to the input signal. Essentially, a measurement of how well the output can replicate the input.


Detecting a frequency swept sinusoid and its parameters?

You could compute the instantaneous frequency of the signal in the frame. This can be done as outlined in this paper (e.g. Eq.(9)). You need to compute the analytic signal using a Hilbert transformer:

$$x_a[n]=x[n]+jx_h[n]\tag{1}$$

where $x_h[n]$ is the Hilbert transform of $x[n]$. The instantaneous frequency is given by the (discrete-time) derivative of the phase of (1). If you use a first order difference for approximating the discrete-time derivative you'll get the formula given in Eq. (9) of the paper cited above. Given an estimate of the instantaneous frequency, you could simply fit a line to estimate the sweep parameter.

This little Matlab/Octave script shows how this approach would work with a toy problem:

% create chirp signal
N = 256;
n=(0:N-1)';
fx=(0.1+0.01*n/N);
x=sin(2*pi*fx.*n);

% compute analytic signal
xh=hilbert(x);
xr=real(xh);xi=imag(xh);

% compute instanteneous frequency (could be done over fewer points)
tmp1 = xr(1:N-1).*xi(2:N) - xr(2:N).*xi(1:N-1);
tmp2 = xr(1:N-1).*xr(2:N) + xi(1:N-1).*xi(2:N);
f=atan2(tmp1,tmp2)/(2*pi);  % (biased) estimate of instanteneous frequency

% linear regression over center values of f
skip = floor(N/4);    % skip unreliable points at beginning and end of f
na = skip:(N-skip);
A=[ones(length(na),1),na(:)];
u=A\f(na);

f0=u(1); k=u(2)/2;  % remove bias
[f0,k]

% reconstructed linear frequency
fr = f0 + k*n;

plot(n,fx,(2:N),f-k*(2:N)',n,fr,':')
legend('exact','estimated','reconstructed')
axis([0,N,f0*.9,(f0+k*N)*1.1])

The plot shows the actual linear frequency of the signal (blue), the estimate from the analytic signal (green), and the reconstructed frequency obtained from fitting a line through the estimated frequency. Note that the estimate is unreliable at both ends, so the line is fitted to the center values of the frequency estimate.


How to design bandpass filter in python when centre frequency is greater than sampling frequency?

As far as the digital signal is concerned, the signal only exists uniquely in one Nyquist interval (which is $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate. Every other interval, such as in your range of interest that is well above the sampling rate, is a replica of the primary interval.  Given that, note that the range from $f_s/2$ to $f_s$ is the same as $-f_s/2$ to $0$, so we can also describe the primary Nyquist interval as being from $0$ to $f_s$ (and Matlab's fftshift command converts a frequency axis from $0$ to $f_s$ to $-f_s/2$ to $+f_s/2$ according to this relationship). The following picture helps illustrate this:

This is formed because sampling is multiplying the time domain signal by an impulse train. An impulse train in time is an impulse train in frequency (as represented by the second row in the figure above). Multiplying in time is convolving in frequency, so the result is to convolve our analog spectrum with the spectrum of the sampling process (an impulse train) resulting in the replicas of the spectrum in the Digital Output Spectrum as shown. Since digitally every higher frequency is an exact replica of what is given from $-f_s/2$ to $+f_s/2$, we need only concern ourselves with the frequency over that range, which is the unique freqency span of a digital signal.  However what I refer to as the "unwrapped view" of the digital signal has helped give me a much more intuitive insight into many multirate signal processing challenges, as well as undersampling (or bandpass sampling) as discussed here. The "unwrapped view" is done by extending the digital frequency axis from $-\infty$ to $+\infty$ while maintaining mathematical equivalence, which means replicating the same signal in each unique span of $n[0:f_s]$ where n is an integer extending to $-\infty$ to $+\infty$ as I have done in the figures here for the "Digital Output Spectrum".

That said, you would need to map the signal of interest to the primary interval and then design the filter in that space.

Further if your signal is real, the positive and negative frequencies would be complex symmetric (as diagrammed also in my figure), in which case we can just concern ourselves with the signal from $0$ to $f_s/2$. I suspect but have not confirmed the range in the python tool you referenced, given as 0 to 1 refers to DC to $f_s/2$ specifically.

As you suggest, we can do bandpass sampling as long as the A/D Converter has an analog input bandwidth that supports this, and I won't get into the detail of it here except to say that the jitter requirements on your sampling clock are more stringent with bandpass sampling. This implies you have an appropriate bandpass filter (refered to as the anti-alias filter) in your analog system prior to your A/D Converter, as the A/D Converter will downconvert every signal in every Nyquist zone if unfiltered energy exists (and within the analog bandwidth of the A/D Converter, which effectively becomes part of the anti-alias filter). This bandpass filtering operation is depicted in the figures below, showing how the energy in the other Nyquist zones gets translated to the same digital spectrum. Without analog bandpass filtering, all these bands will downconvert through the sampling process to the same frequency range, leading to a significant increase (increase in the noise floor even if there are no other signals present if the amplified noise floor in front of the A/D Converter is not bandpass or lowpass filtered in any way, as said earlier, out to the input analog bandwidth of the A/D Converter.

That said, in your case of a 100 MHz center frequency, the 100 MHz lands directly on a harmonic of the sampling clock and therefore maps to 0. Unless you are dealing with a complex signal or a real signal in a synchronous system, this would be a poor choice for a center frequency or sampling rate (by synchronous, I mean the signal of interest is coherent to the sampling clock such that you are assured it is exactly centered on 100 MHz and complex symmetric on each side of 100 MHz). The reason is, any offset in center frequency will cause aliasing in the image that gets mapped to 0 resulting in severe distortion. (Exceptions to this would be in closed-loop systems that will drive any offset to exactly 0). If you use complex sampling with resulting quadrature I and Q signal paths then the spectrum is unique from $-f_s/2$ to $+f_s/2$ and this will be a non-issue. For real signals it is common practice to center the signal near $f_s/4$ and multiples there-of just as I have done in these figures as it optimally relaxes the requirement on your analog anti-alias filter (the images will be equidistant as in my figures). It is common to then use a complex down-converter to then move the signal to baseband appropriately without aliasing issues. So to do this you would either shift your center frequency to 102.5 MHz, or shift your sampling clock (using some math) such that 100 MHz is mid-way in a Nyquist zone. You would then design a bandpass digital fitler that is in your primary Nyquist zone of $0$ to $f_s/2$ or as your tool is givng you in normalized terms $0$ to $1$.


From symbol rate to sample rate fractional ratio

The symbol rate is independent of the sampling rate. You can sample the symbols at any desired rate that is greater than twice the bandwidth of the symbol, where higher rates would simplify the required analog filtering after the DAC (so higher sampling is better in this regard and the initial ratios the OP is using are reasonable).

The bandwidth of each symbol is driven by the pulse shaping that may be used. With typical pulse shaping such as Raised Cosine, two samples per symbol would be a reasonable minimum sampling rate to use but would also require tighter analog filtering.

To see this, consider each “symbol” as it would appear as an analog waveform and observe how you can sample that any any rate as long as you observe Nyquist’s criteria to sample greater than twice the highest bandwidth of interest in the signal. Going from digital to analog with a DAC is this same process in reverse, in that how many samples do you need to properly reconstruct the desired waveform.

Assuming this design starts with symbols sampled at one sample per symbol at a 3 MHz symbol rate, typically this would be interpolated (upsampled by inserting zeroes and then filtered) where the filter would be the pulse shaping filter defined by the waveform transmitter specifications (often a root-raised cosine filter but must be defined).  This upsampling could be matched to the transmitter sampling DAC frequency using polyphase resampling or done in stages.  If in stages resampling to only 4 MBps would not be sufficient for pulse shaping given the explanation above, but it would be reasonable to resample to 10 MHz for pulse sampling and followed by 2x interpolation for the final DAC sampling.


Cepstrum: IFFT or FFT to enter the cepstral domain?

For Cepstrum I always have used to this steps:

Apply hamming windows in the signal
Apply FFT
Get magnitude
Convert to log scale
Apply IFFT

The equation for cepstrum:

IFFT(log(abs(FFT(s))))

But you can use FFT or IFFT, take a look:

IFFT(log(abs(FFT(s)))) == real(FFT(log(abs(FFT(s)))))

The difference is the scale representation, if do you end using FFT you need extract just the real information, for both above equations you will get the same shape:

For IFFT(log(abs(FFT(s)))):

For real(FFT(log(abs(FFT(s))))):

This is a cepstrum shape example from 4096 points sine in 440hz sampled at 44100hz


Are oversampling and downsampling the same thing?

OK. it looks like I got things mixed up.

The book on DSP I have mentions oversampling and that it is / can be followed by downsampling/decimation, ... and therefor I (wrongly) assumed that "oversampling" included both steps.

Mystery solved.
Thanks for the correction.

Kristoff


Are the frequencies a histogram in the frequency domain?

No.

A histogram is the visualisation of a sequence of discrete counts. For example, a population pyramid can tell you the distribution of ages in a population.

The frequency spectrum represents strength. It doesn't have to be integer and it does not really mean that it represents counts of sinusoids at a specific frequency. The most trivial counter example of this is what is known as the Direct Current (DC) value which is frequency 0 (Zero) Hertz.

Although, you could possibly formulate a view of the frequency spectrum as some form of histogram (counts of elementary parcels of sinusoids), it would be very difficult to fit this view to the average value across the sample (which is what 0 Hz represents).

So, the frequency spectrum is not a histogram, it has nothing to do with a histogram. But what is it?

It is difficult to tell you without going through equations and properties because they represent our understanding of how the world works (for example, linearity and how things can be decomposed to parts and recomposed from parts). But here is an attempt:

This is a string trio:

Notice that the instruments have more or less the same shape but their size changes. Now, without going into equations, take it as a rule of thumb: The bigger the instrument the "lower" (in frequency) is its sound. (And, they can get...big).

When they all play together, they make up one continuous stream of sound that hits one's ears or a microphone.

The accurate analogy here is that the frequency spectrum of the sound recording would tell you how loud was each one of these 'violins' playing.

Joseph Fourier worked out a finer analogy that is very useful in DSP and goes by the name Discrete Fourier Transform (DFT).

He postulated that given any periodic waveform, there is a way to decompose it down to an infinite sum of sinusoids of varying strength and phase.

This is what the DFT does and from the DFT we derive various frequency spectra (namely, amplitude, power, phase).

  ...since if there are two signals of the same frequency but different amplitude, wouldn't they have a representation on the frequency graph where it's not possible to distinguish between them?

You see how this histogram view of the frequency spectrum doesn't really work now I hope. Two signals of the same frequency contribute to the strength of the same frequency in "the mix" and therefore do not appear as two individual measurements.

  I'm even more confused because this example shows 20 millisecond samples being taken from a Fourier spectrum...

But what you show is not a frequency spectrum. What you show is a spectrogram.

Without going into equations: The DFT decomposes a periodic waveform (any waveform) into an infinite sum of sinusoids. Because of the properties of those sinusoids, the DFT can tell you that a sinusoid is present in the mix but it cannot tell you WHEN it was present.

So, if you had a 1 hour recording of just two violins, a big one and a small one but with the small one just performing for 2 seconds, the DFT will still show high frequencies but it will not be able to inform you WHICH 2 SECOND SEGMENT of the 1 hour recording was it that the small violin was louder.

For this reason, the Short Time Fourier Transform (STFT) was developed. What this does is decompose the signal down to small (possibly) overlapping windows and obtain the DFT's of those windows.

This leads to an intermediate representation of the signal where we now can tell (within limits) which frequencies were present at which times.

This is what the last diagram in your question is depicting and this is why it still retains an axis of "Time". In that diagram, Time runs left to right and frequency runs bottom to top.

Hope this helps.

PS: I think that I understand why you don't feel like going through the equations but if you make a start, first with understanding the mathematical notation and then towards reading and understanding what it is that the equations are trying to describe, then I think that you will be greatly benefited in the long term.


The simplest approach would be to divide the gradient rather than subtract it. Here's what the result looks like:


In traditional FSK, the frequency separation between pulses is $1/T$, where $T$ is the symbol rate. This guarantees orthogonal pulses while maintaining a continuous phase. The phase is continuous because a pulse with the next higher frequency traverses exactly one more cycle than the one with the lower frequency (during one symbol interval).

However, you can easily show that in order to guarantee orthogonality between pulses, a frequency separation of $1/2T$ is actually sufficient. This is indeed the minimum frequency separation guaranteeing orthogonality. Now a higher frequency pulse traverses only half a cycle more than the pulse with the next lower frequency. The smaller frequency separation reduces the bandwidth compared to traditional FSK, whereas the probability of error remains unchanged due to the orthogonality of the pulses. This is why MSK is often referred to as bandwidth-efficient (compared to traditional FSK).

Note that for MSK to have a continuous phase the initial phase of a pulse must be either $0$ or $\pi$, depending on the phase of the previous pulse (because some pulses have a phase value of $\pi$ at the end of a symbol interval due to the minimum frequency separation of half a cycle per symbol interval).

This additional "phase encoding" inherent in MSK can actually be used in the receiver by applying the Viterbi algorithm. In this way, MSK can perform better than traditional FSK. It can overcome the 3dB loss of FSK as compared to PSK.

Also note that MSK can be viewed and implemented as a special case of offset keyed PSK (with sinusoidal pulse shapes).


Difference between acoustic emissions and ultrasonic waves

The propagation of a mechanical wave within a medium depends on the characteristics of that medium (e.g. density).

An easy way to experience this is to knock a piece of material with a sharp hit and listen to it trying to dissipate the energy you gave it. Consider for example a bell. When you hit it, it goes bonnnnnnnnnnng at some tone. If you were to weld a few rings around the bell and try to hit it again, it would go bonng. Different density / stifness of the material, different sound.

Acoustic emissions are akin to the sound of the bell and the cause of the "ringing" are the internal cracks. Consider for example the sound of ice cracking (similar sounds are produced even as ice cracks on its own, without forcing it externally) or even closer to your question, glass cracking under thermal shock.

Within a block of ice, the acoustic emissions are low enough in frequency that can be heard. But, if that was a piece of kevlar breaking, the disturbance would propagate with different characteristics through the material. You can see that we are already at hundreds of kHz in this case and we would be needing special equipment to sense the material vibrating.

Hope this helps.


The answer to your question "is it a problem" is really up to you to determine. Sometimes the simple things work well enough.

There is a technique called Dynamic Time Warping:

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

and Matlab has a function dtw

https://www.mathworks.com/help/signal/ref/dtw.html

which can use to figure out if it is a problem.

I suggest you pick a reference interval, that would be the closest to the average shaft interval and use that to warp the rest against.


Techniques to reduce delay in detecting face from video using vision.CascadeObjectDetector in matlab

There are a few of things you can try:

Definitely move FaceDetect = vision.CascadeObjectDetector; outside of
the loop. You only need to create the face detector object once.
Re-creating it for every frame is definitely your performance
bottleneck.
vision.VideoFileReader returns a frame of class 'single' by default.
If you change the output data type to 'uint8', that    should speed
up    the face detector.Use obj=vision.VideoFileReader('basu.avi',    'VideoOutputDataType',
'uint8');
vision.VideoFileReader can also do the conversion to grayscale for
you. Use obj=vision.VideoFileReader('basu.avi',
'VideoOutputDataType', 'uint8', 'ImageColorSpace', 'Intensity'); This
may be faster than calling rgb2gray.
Try limiting the size of the faces being detected using 'MinSize' and
'MaxSize' options of vision.CascadeObjectDetector    and/or try
downsampling the frame before detecting faces.


What should I do know to know the frequency from complex signal?

When you say "complex" signal, I assume you mean just a signal with many frequencies, and not a "complex-valued" signal, in which the the signal values themselves are complex numbers. Let me know if I am incorrect in this assumption. (Actually, computation of DFT is the same in either case).

I'm also assuming your signal is discrete.

The discrete Fourier Transform, $F(\omega)$, can be found by solving the following equation (where N is the number of samples your signal, $f[n]$, contains)

$$F(\omega) = \sum_{n=0}^{N-1} f[n] e^{\frac{-j 2\pi \omega n}{N}}$$

Alternatively, you could use the function fft() on MATLAB.


I do not understand this simple bandstop filter example

fir1 designs your filter using a Hamming window one bin larger than the degree you specify. 41 bins is too small to give you the frequency resolution you want. Try a much larger degree, say 300-600.


How to get the autocorrelation of this finite discrete function?

Calculating the autocorrelation of this function is relatively easy. You just need to follow the formula you have written:

\begin{equation}
\ r_{xx}(l) = \displaystyle\sum_{n=-\infty}^{\infty} x(n)x(n-l) \end{equation}

It gets much easier because your function is $\neq0$ for only a few values ($2N+1$ exactly).

Because of this the previous sum can be replaced with this much simpler one:

\begin{equation}
\ r_{xx}(l) = \displaystyle\sum_{n=-N}^{N} x(n-l) \end{equation}

Which is basically just counting how many times $x(n-l)$ is not zero on that finite interval. Because of this it appears easily that

\begin{equation}
r_{xx}(l)=max(2N-1-|l|,0)  \end{equation}

You almost had it in your question but the sign was wrong.


Problems such as those are closely related to "inverse problems", whereby one attempts to reconstruct a signal from an under-determined system of equations, (ill-posed problems) given some apriori information. As @Hilmar mentioned, you cannot reconstruct phase information from the absolute magnitude of a DFT, unless you have prior information that constrains you. For example, in the paper you mention, the author states:

  Both the problem of phase
  retrieval from two intensity measurements (in electron microscopy or
  wave front sensing) and the problem of phase retrieval from a single
  intensity measurement plus a non-negativity constraint (in astronomy)
  are considered, with emphasis on the latter. It is shown [...]

The bold is mine. Even in the paper you allude to, there are constraints that then allow for a Maximum Likelihood framework to be applied, or one where sparsity is rewarded, etc.

  Edit: The prior information I have is : Time-domain signal is
  positive-valued. Time domain signal is zero for in some points, i.e.
  signal has a support P⊂{0,1,2,⋯,N−1}

All your signal being positive valued is irrelevant to the phase information of the entire DFT spectrum. To see this, consider this following:

Note that having an all positive signal $x[n]$ means you are guaranteed to have a DC bias. This means that the first bin in the absolute magnitude of you DFT corresponding to $|X(0)|$ will be strictly positive. Let us split the input signal into its zero-mean partition, ($z[n]$) plus the mean value, $m$. Then:

$$x[n] = z[n] + m$$

The DFT of the result is, (exploiting the linearity property of the DFT):

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

Recall that the DFT is nothing but a projection of your signal $x[n]$ onto an orthogonal basis set composed of complex exponentials. Save for the DC component at $k=0$, all the complex sinusoids have mean $0$. Thus the second term $\sum_{n=0}^{N-1} m \ e^{-j \frac{ \Large 2 \pi n k }{N}}$ completely vanishes for all $k \ne 0$. Hence, we can now write the DFT of $x[n]$ as being:

$$
X[k] = \left\{
  \begin{array}{l l}
     \sum_{n=0}^{N-1} z[n] \ e^{-j \frac{ \Large 2 \pi n k }{N}} + mN = mN & ,k = 0\\
    \sum_{n=0}^{N-1} z[n] \ e^{-j \frac{ \Large 2 \pi n k }{N}} & , k \ne 0
  \end{array} \right.
$$

This then shows that your signal always being positive is completely irrelevant to the phase signal spectrum.

Hence, this prior information ($x[n] \ge 0, n \in 0, 1, ... N-1$) is irrelevant, and no recovery is possible given $|X[k]|$, with only this constraint.

You are right in that these types of problems have some solutions to them, but it is also important for you to also remember that you have to define your constraints before you can solve for your unknown entities. That is, you have to give prior information that you can utilize.

As a simple example, given the absolute magnitude vector $|X(f)|$, there is an infinite number of possible solutions to what values $\angle X(f)$ can take on. However, if you "know" before hand that, say, all the imaginary co-efficients of $X(f) \in \mathfrak{C} $ are zero, then we can easily recover what the real co-efficients are up to a sign ambiguity.

The take home message however, is what prior information do you have, and how is it related to the DFT. Hope this helped.


What is the relationship between poles and system stability?

I agree with Peter K.'s answer, but I would like to add one important point: the two statements in the question are only true for causal systems. The most general statement about stability for LTI systems described by rational transfer function is:

  An LTI system with a rational transfer function is stable if the region of convergence (ROC) of its transfer function includes the $j\omega$-axis (for continuous-time systems), or the unit circle $|z|=1$ (for discrete-time systems).

Since for causal systems the ROC is to the right of the pole with the largest magnitude (or outside the circle with radius equal to the largest pole radius in the $z$-plane), the above statement coincides with the statements in the question. However, non-causal systems can have poles in the right half-plane (or outside the unit circle of the $z$-plane) and still be stable.

Take as an example the discrete-time transfer function

$$H(z)=\frac{z^2}{(z-a)(z-b)},\quad a<1<b\tag{1}$$

which has one pole inside the unit circle and one pole outside the unit circle. You can't say that the corresponding system is unstable, because it depends on which system you mean. The transfer function (1) has 3 possible ROCs, corresponding to 3 different systems:

$|z|>b$
$|z|<a$
$a<|z|<b$

The first one corresponds to a causal system which is unstable, because one of its poles is outside the unit circle (the ROC does not include the unit circle). The second one corresponds to an anti-causal system which is also unstable, because it has one pole inside the unit circle (again, the ROC does not include the unit circle). Finally, the third ROC includes the unit circle, i.e. it corresponds to a stable system. It is a non-causal system with a two-sided impulse response which decays exponentially as the time index $n$ goes to $\pm \infty$.


Why downsample exhibit wider bandwidth when doing discrete time fourier transform?

You don't get "more" bandwidth. That's simply a misunderstanding of your x-axis. Your signal after downsampling occupies a $M$ times larger part of the Nyquist rate, but that's because the Nyquist rate got reduced by a factor of $M$, not because the signal got wider!


Quantization is a process of representing real values of a signal using integer numbers. That process results in loss of accuracy, which is determined by the selection of the smallest value of quantization level.
If the smallest quantization level is $\Delta$, then every signal level will be rounded to the multiple of that $\Delta$. For example, if a signal level is $x=3,3 V$ and $\Delta = 1V$, then the corresponding quantized representation of the signal level will be $Q(x)=3\Delta$ and the resulting $\hat x = 3_{10} = 11_2$. If $\Delta = 0.5V$, then the signal level will be $Q(x)=7\Delta$ and the resulting $\hat x = 7_{10} = 111_2$.
So, to rephrase the original point: the weight of the least significant bit of the binary code word equals $\Delta$.


Typically "smoothing" means "replace the current value with average over the neighboring ones". Most common is energy smoothing, where the smoothing results in the energy average over the smoothing interval and the phase information is lost. Complex smoothing can be done as well but it's tricky business because of phase wrapping.

Energy smoothing can be expressed as
$$Y(k)=\sqrt{\frac{1}{N}\cdot \sum_{i=0}^{N-1}X(i)\cdot X^{*}(i)\cdot W_{k}(i)}$$

where $W_{k}(i)$ is some suitable window function. In the case of, say, third-octave smoothing this could be derived as the the magnitude squared of the transfer function of a third octave band pass filter around frequency k. This also means that for a, say, 1024 point FFT you need to design 1024 different bandpass filters, so that's a fair bit of work.

Things can be simplified if the exact shape of the smoothing filter is flexible. Rectangular smoothing can be done as
$$Y(k)=\sqrt{\frac{1}{b-a+1}\cdot \sum_{i=a}^{b}X(i)\cdot X^{*}(i)}$$

where
$$a = round(k*2^{-\frac{1}{2\cdot n}}), b = round(k*2^{\frac{1}{2\cdot n}})$$

are simply the indices of the band edges for $n^{th}$ octave smoothing.

There are a few more methods that are in between in the arbitrary window and the rectangular one in terms of complexity.


Apply Kalman filter to remove process noise, given zero measurement error

You're misunderstanding how the Kalman filter framework works. Your measurements are a reflection of what the underlying state of the system is at each time instant. If you have zero measurement covariance, then your measurements perfectly reflect the system state, so there's no need for any filter whatsoever. It's not likely that this is truly the case.

Process noise is a representation of uncertainty in the model for how your system's state changes as time goes on. It's not something to be "removed" by the Kalman filter. Instead, it gives you the ability to accommodate higher-order effects that you aren't modeling in your state transition matrix.


How do the motion vectors work in predictive coding for MPEG?

To simplify your confusion - there two processes :

1. Motion estimation
2. Motion compensation

Before we talk about the estimation, we should talk about the Motion compensation.

Let say, the $ Image_{t}(x,y) $ is split in block $ Blocks_{t}[k](x',y') $.

The task of Motion compensation is to produce $ Blocks_{t}[k](x',y') $ from any region of
$ Image_{t-1}(x,y) $.

Hence another block not necessarily aligned at 16x16 boundary is a best possible match $ Blocks_{t-1}[k](x'+mx,y'+my) $

Here, $ mx, my $ is called motion vectors.

We can calculate the error between the target and reference as

$$ Err_{t}[k](x,y) = Blocks_{t}[k](x',y') - Blocks_{t-1}[k](x'+mx,y'+my) $$

So now, the encoder basically transmits $Err_{t}[k](x,y)$ (with DCT and quantization) and $ {( mx, my) }[k] $  for each block,.

So encoder has 2 work to do:

1. Motion Estimation
The process or estimating $ { mx, my }[k] $ for every $ k $ such that $ Err_{t}[k](x,y) $ is minimized is called Motion estimation.

2. Generation of error image after Motion Compensation
The process of constructing $ Blocks_{t}[k](x',y') $ from $I_{t}$ image pixels and $ { (mx, my) }[k] $ is called Motion compensation.
The error image is what gets transmitted.

Finally, decoder can re-do motion compensation on its' own using moiton vectors and the error image to make final re-construction of the image.

Now we realize a few points:

Best Motion estimation helps minimize the energy that is required to be transmitted and hence optimizes the bits for a given quality.
However, even if $ { (mx, my) }[k] $ is not ideal or if the scene has significant change over last picture, the $ Err_{t}[k](x,y) $ is always transmitted to receiver - hence the reconstruction is always perfect (modulo the loss created by the quantization). Hence, even if you have suboptimal motion vector or redundancy is not much, the reconstruction is still always perfect albeit with more bits!
Every block $ Blocks_{t}[k](x',y') $ is motion compensated on its own right - hence even if the actual motion vectors of any neighboring blocks has no effect in the construction. Hence, it is not necessary to have motion vectors perfectly aligned for perfect reconstruction to be possible.
Albeit algorithm  exists who are smart enough to guess that if $Blocks_{t}[k]$ has motion vector $( mx ,my )[k]$ the guess for $Blocks_{t}[k+1]$ might just be closer to that only.
Finally suppose the next picture is completely different it is possible that energy for $$ Energy (Err_{t}[k](x,y))  > Energy ( Blocks_{t}[k](x',y') ) $$.

In such cases, it might be more advisable to transmit block directly without prediction than to send the difference. This is also possible in the encoder by a provision called INTRA block.


Multiplier in digital signal processing?

Once upon a time, chips (processor ICs) did not have anywhere close to millions (now billions) of transistors on them.

A typical fast multiplier, say made out of a tree of layers of carry-save-adder cells, requires a lot more logic gates (thus transistors) than a simple adder, thus also requiring a lot more power, and with less of them fitting on a particular silicon die area.

Repeated conditional-add-and-shift is slower (in the same transistor technology) but requires a lot less logic gates (thus less transistors) when done sequentially.

If an (ancient) FPGA has tons of small memories (but no multipliers), the memories might be useful as look-up tables that are faster than the equivalent logic equations implemented in fabric connected LUT cells.


Confusion with regards to units in Sinc pulse

The length of the pulse, $T$, would be in units of time (for example, in seconds). The bandwidth would be in units of frequency (for example, in Hz). So

$$[T][\Delta f]=sec\frac{1}{sec}=1$$

That means that $Z$ is a dimensionless quantity. That makes sense, as $Z$ is just a natural number that describes the number of times that something happens (in this case, the number of zero crossings), and should not have units at all.


Calculating smoothed derivative of a signal by using difference with larger step=convolving with rectangular window

This is a hard question to handle generally. Smoothing with a rectangular window is used all the time (often called a "moving average"), so that's not necessarily a problem. I'm not sure what ringing you're referring to, perhaps the sidelobes of the rectangular window's frequency response.

Differentiation is inherently a highpass operation; the ideal continuous-time differentiator has a transfer function of:

$$
H(s) = s
$$

So that its magnitude response is:

$$
|H(j\omega)| = \omega
$$

The gain of a differentiator therefore increases monotonically with frequency. If your signal contains high-frequency noise, then it could be amplified by applying a differentiator. To combat this, two approaches are apparent:

Design a more sophisticated differentiator filter that has the desired linear magnitude response over the portion of the band that covers your signal of interest, then sharply attenuates higher frequencies. You could design such a filter using least-squares methods or the frequency-sampling method, for example.
Use a cascaded approach where you first suppress all of the high-frequency noise that you can using a lowpass filter, then follow with the differentiator. The frequency coverage of the differentiator need not be as tight since the lowpass filter will eliminate the out-of-band noise.

The methods should be roughly equivalent if you're using linear filters; you could think of the first single-filter approach as just a cascade of a differentiator and a lowpass filter. As you noted, the central-difference approach can be modeled in this way. It's hard for anyone to say that it's "bad" without any knowledge of your application. My main thought would be that it is "bad" if the smoothing operation tangibly attenuates your signal of interest, such that the derivative estimate is no longer useful. However, if the parameters of the signal are such that you can smooth the noise without noticeably distorting the signal (i.e. if the signal is well-oversampled), then it could be a win.


A continuous wavelet frame can be discretized with perfect reconstruction under some technical conditions. If this is possible, you obtain a discrete wavelet frame, often with some amount of  redundancy (depending on the wavelet shape). This is a result formalized by I. Daubechies, for instance.

For a long time, apart from the Haar or the Franklin wavelets, little was known about the possibility of obtaining a critical (non-redundant) discrete wavelet transform. Until the inception of the Meyer wavelet, and the subsequent ones (Daubechies, Symmlets, Coiflets, etc.). Their somewhat weird shapes is the sign that not every wavelet shape is possible. This is even truer as the complex Morlet is complex, hence somewhat already twice redundant for real signals.

So, Morlet cannot be critically discretized with 2-band (dyadic), as you already guessed, as you wanted to compare CWT and DWT. However, multi-band  filter banks are an option to build orthogonal counterparts of "Morlet"-like design. Malvar wavelets, Modulated complex lapped transforms or Lapped orthogonal transforms are  instances (see a Panorama of wavelets for references), and their can be  further decomposed across scales as $M$-band wavelets.

Or you can accept a non perfect reconstruction dyadic decomposition akin to Morlet.


How do I find the Energy Density Function of $g(t)$ if i am not given an input or impulse response?

If you're referring to energy spectral density, then you need to compute the squared magnitude of the Fourier transform of $g(t)$:
$$\big|G(\omega)\big|^2=\left|\int_{-\infty}^{\infty}g(t)e^{-j\omega t}dt\right|^2\tag{1}$$


Your question is very broad and I won't go read the papers and implement it for you. If you could show specific problems, we might help you better (or the guys at stackoverflow for programming questions).

Meanwhile, you can have a look at these two resources, which treat a similar problem, including some code:

http://www.pyimagesearch.com/2014/09/01/build-kick-ass-mobile-document-scanner-just-5-minutes/

https://mzucker.github.io/2016/08/15/page-dewarping.html


Equivalence between noise covariance of real and complex models of received signal

Even though I doubt your noise model, here's how you can do it:

Given the complex noise vector $\overline w$, with the covariance matrix $E[\overline w\overline w^H]=\sigma^2\overline{G}$.

Let's separate $\overline w$ into its real and imaginary part (I use a different notation, because I think this is faster to write and read): $\overline w=w_R+jw_I$, where $w_R, w_I$ denote the real and imaginary parts of $\overline w$. Insert this into the covariance definition:

$$\begin{align}E[\overline w\overline w^H]&=E[(w_R+jw_I)(w_R^T-jw_I^T)]\\&=E[w_Rw_R^T+w_Iw_I^T]+2jE[w_Rw_I^T]=\overline{G}\end{align}$$

Let's compare this to the expression we get, when explicitely separating the real and imaginary parts:
$$\begin{align}E[ww^T]&=\begin{pmatrix}E[w_Rw_R^T] & E[w_Rw_I^T]\\E[w_Iw_R^T] & E[w_Iw_I^T]\end{pmatrix}=X\end{align}.$$

Hence, in order to define $X$ (which is the matrix you are interested in), we need to know $E[w_Rw_R^T], E[w_Rw_I^T]=E[w_Iw_R^T], E[w_Iw_I^T]$.

Let's do the comparison to $E[\overline w \overline w^H]$, which we have calculated before. We know that $G_I=0$ (i.e. the imaginary part is zero). Hence, $E[w_Rw_I^T]=0$, i.e. real and imaginary part are uncorrelated. What we dont know yet are $E[w_Rw_R^T]$ and $E[w_Iw_I^T]$. We only know that $E[w_Rw_R^T]+E[w_Iw_I^T]=\overline G$.

How to solve it? There are two possibilities:

1) (The more common assumption) We assume that real and imaginary parts of the noise have equal variance. This is often reasonble, e.g. when the noise is baseband noise and was obtained from downconversion a passband signal. Then, we have $E[w_Rw_R^T]=E[w_Iw_I^T]=\frac{1}{2}\overline{G}$.

2) Real and imaginary parts are not i.i.d. Then, the noise is not proper and the pseudocorrelation $T$ of the noise does not vanish. It is defined by $T=E[\overline w \overline w^T]$ and we calculate it to be (assuming $E[w_Iw_R^T]=0$)
$$T=E[\overline w \overline w^T]=E[w_Rw_R^T]-E[w_Iw_I^T]$$
In order to accurately describe a non-proper noise, the value of $T$ must be known. Then, we have two unknowns ($E[w_Rw_R^T]$ and $E[w_Iw_I^T]$) and two equations $$\begin{align}\overline{G}&=E[w_Rw_R^T]+E[w_Iw_I^T]\\T&=E[w_Rw_R^T]-E[w_Iw_I^T]\end{align}$$
from which you can get the unknown via $E[w_Rw_R^T]=\frac{1}{2}(\overline{G}+T)$ and $E[w_Iw_I^T]=\frac{1}{2}(\overline G-T)$.

Given this derivation, I still doubt your noise model. To me it does not make any sense, that in the equation $y=\overline G \overline a+\overline{w}$, the noise should have variance $\sigma^2\overline{G}$ due to ISI. Can you ensure that $\overline G$ is Hermitian?


strange phenomenon in spectrogram of Windowed Fourier Transform [edited]

This happens because your window is too short. I don't have access to a plotting tool right now, but imagine for a second a slowly varying sinusoid that you chop up into pieces, and these pieces are shorter than the period of your sinusoid. If you take the Fourier transform of each of these pieces, some of your chunks will capture more energy of this sinusoid than others because they will contain crest or through of the function, while others will contain regions close to 0.

If you increase your window size, you'll be able to catch several cycles of your slow sinusoids inside each window, and the periodicity will go away. The tradeoff is that your windows will capture larger lengths of time, so boundaries between events will be more blurred.

If this still doesn't make sense, I'll try to include some plots later.


Assumptions for Hurst exponent calculation

The best is to calculate Hurst for stationary data (as then you can somehow find out if there are long range dependency signals which are not part of AR or MAs) but you can do non stationary as well. You can calculated froma fractal or wavelet transforms. Results are different

In calculating Hurst you have to divide the data in short windows and average the results, (place a range the last 8 elements). Hurst can fluctuate significantly depending on the period you select.