Light Detection and Ranging (LIDAR) Range Bias Compensation with Light Dependent Resistor (LDR)

Ideally, you would like to explain the LIDAR signal as a combination of one or more other variables, here LDR and/or Temperature and then remove the effect of those variables from the LIDAR signal.

The first part can be achieved with linear (hopefully) regression and specifically Multiple Linear Regression where you would be fitting a model like LIDAR = Constant + LDR_Coef * LDR + Temp_Coef * Temp.

With the availability of just the LDR data (code in GNU Octave but easily transferable to other platforms too):

%Load the data
D = dlmread("myfile.txt",";"); %myfile.txt is the file you provided which is delimited by semicolon.
% Adjust the LDR signal
D(:,3) = log(D(:,3));
% Fit the model using regress.
% For more information see: https://octave.sourceforge.io/statistics/function/regress.html
[B, BINT, R, RINT, STATS] = regress(D(:,2), [ones(size(D(:,3))), D(:,3)]);

After you do this, you get your model coefficients in B and a few other things (more about this later on).

It is now possible to express LIDAR as a function of LDR:

% Plot LIDAR and LIDAR as a function of LDR
r=0:(length(D)-1);
plot(r,D(:,2),B(1)+B(2).*D(:,3));
xlabel('Time');ylabel('LIDAR / LDR');legend('LIDAR','LDR');grid on;

This now looks like:

This concludes with the fitting part. Now we want to remove the variation due to LDR, which is the "easy" part:

LIDAR_adjusted = D(:,2) - B(2).*D(:,3);
plot(r,D(:,2),LIDAR_adjusted);
xlabel('Time');ylabel('LIDAR');legend('LIDAR','LIDAR_adjusted');grid on;

Which results in:

The STATS vector contains a number of coefficients that inform you about how good is this fit. In this case 7.6591e-01   1.4952e+03   0.0000e+00   7.8229e-04. So, the $R^2$ is not ideal but not too bad too (that is the first number, it can be between $0$ and $1$ and it is ideal at $1$), the rest of the numbers are not bad too (the error is very small), but if you look at the residuals (that is the difference between the LIDAR and the LIDAR as a function of LDR) it is not exactly noise. Here is the residual (R):

plot(R);grid on;

Which results in:

There are clearly trends in the residual which means that LDR on its own cannot explain the variation in the LIDAR signal. You can add 'r' as a regressor in there, which would return three coefficients for the model LIDAR = B(1) + B(2)*D(:,3) + B(3)*r and that would "take away" that linear upwards trend as a function of time and improve the error but you still get those two humps because the sensor responses are non-linear.

To conclude with this, if you wanted to add the temperature to the model too, just change the regress line to:

[B, BINT, R, RINT, STATS] = regress(D(:,2), [ones(size(D(:,3)), D(:,3), D(:,4)]);

Assuming that you will add another column for the temperature. The rest of the process remains the same (with reasonable adjustments for removing the components further below).

The other thing that you can do of course is a Discrete Fourier Transform (just on the LIDAR) and retain the first few harmonics from it for detrending. The latter assumes that there is a periodic well defined effect of light (and / or temperature and other parameters) on LIDAR which can be approximated with a set of periodic components. And in that case, the adjustment comes as a function of the time of day "for free".

Hope this helps.


How to filter my noisy square wave so that i can recover my signal using the Chebyshev Filter

Bandpass filtering won't do what you're after.  If I take a noiseless square wave with the spectrum in red, and design a bandpass filter for it (response in blue):

and then apply this filter to the original (light blue) square wave, then I get the dark blue response.

I suspect that a better approach might be to use a median filter.

If I apply some gaussian noise to the square wave (blue line), and then apply a median filter (orange plot) then the result is below.

These plots were generated in python; I don't have easy access to Matlab. Code below.

from scipy import signal
import math
import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0,20*math.pi,1000)
sq = scipy.signal.square(t)

b, a = signal.cheby1(4, 0.1, [0.01, 0.2], 'bandpass')

plt.rcParams['figure.figsize'] = [20, 10]
plt.figure(1)
w, h = signal.freqz(b, a)
w2, h2 = signal.freqz(sq, 1)
plt.plot(w, abs(h)*600)
plt.plot(w2, abs(h2), color='red')

plt.figure(2)
plt.rcParams['figure.figsize'] = [20, 10]
plt.plot(t,sq)
sq_f = scipy.signal.lfilter(b,a,sq)
plt.plot(t,sq_f,color='blue')

sq_noisy = sq + np.random.normal(0,1,1000)
plt.figure(3)
plt.plot(t,sq_noisy)
sq_noisy_f = scipy.signal.medfilt(sq_noisy,21)
plt.plot(t,sq_noisy_f,color='orange')


Non-zero DFT components where zero is expected?

These are just rounding errors, and you only get a flipped sign for phase values of $\pm\pi$, which is not an issue because the phase value is ambiguous, i.e., adding or subtracting multiples of $2\pi$ doesn't change the value of the complex coefficient.
I assume that your code is just for an educational purpose, because it is neither efficient nor numerically accurate.


When convolving two functions that are constants in a region and 0 everywhere else, where does the integration start?

For simplicity, let's talk about one-dimensional convolutions only.

You need to understand that the result of a convolution is a function and the alleged convolution integral is actually an infinitude of integrals, one for each time instant.  By this I mean that to compute the output $y(t)$ of a system with impulse response
$h(t)$ to an input signal at time $t = 5$, we write
$$\begin{align*}
y(5) &= \int_{-\infty}^\infty h(\tau)x(5-\tau)\mathrm d\tau,\\
\end{align*}$$
and compute $y(5)$.  For $t = 6.5$, we write
$$\begin{align*}
y(6.5) &= \int_{-\infty}^\infty h(\tau)x(6.5-\tau)\mathrm d\tau,\\
\end{align*}$$
and compute $y(6)$.  For $t = 10.2$, we write
$$\begin{align*}
y(10.2) &= \int_{-\infty}^\infty h(\tau)x(10.2-\tau)\mathrm d\tau,\\
\end{align*}$$
and compute $y(10.2)$.and so on.  After doing this a while, some smart-ass
said, "Hey Ma, I am beginning to see a pattern here!  Lookit, I can
write
$$\begin{align*}
y(t) &= \int_{-\infty}^\infty h(\tau)x(t-\tau)\mathrm d\tau, ~\text{for all} ~t\\
\end{align*}$$
and the world was never the same again.

Suppose that $h(t)$ is $0$ for $t < 0$ or $t > 12$.  Then all the integrals
exhibited above simplify in the sense that the limits can be changed to
$0$ and $12$ since the integrand, being a multiple of $h(\tau)$, is $0$
when $\tau < 0$ or $\tau > 12$.  Suppose in addition that $x(t)$ is nonzero
only when $|t| < 3$.  Now look at the first integral (for $y(5)$ with the
limits changed to $0$ and $12$).  At $\tau = 0$, the integrand is
$h(0)x(5-0) = 0$ since $x(5) = 0$.  As $\tau$ increases towards $12$, the
argument of $x$ in the integrand decreases from $5$ downwards till the
argument hits $-7$ at $\tau = 12$.  Thus, the integrand is nonzero
only for $\tau \in (2,8)$, having value $h(2)x(5-2) = h(2)x(3)$
at one end of the interval and value $h(8)x(5-8) = h(8)x(-3)$ at
the other.  We conclude that
$$y(5) = \int_{-\infty}^\infty h(\tau)x(5-\tau)\mathrm d\tau
= \int_{0}^{12} h(\tau)x(5-\tau)\mathrm d\tau
= \int_{2}^{8} h(\tau)x(5-\tau)\mathrm d\tau.$$
Similarly,
$$y(6.5) = \int_{-\infty}^\infty h(\tau)x(6.5-\tau)\mathrm d\tau
= \int_{0}^{12} h(\tau)x(6.5-\tau)\mathrm d\tau
= \int_{3.5}^{9.5} h(\tau)x(6.5-\tau)\mathrm d\tau.$$
However,
$$y(10.2) = \int_{-\infty}^\infty h(\tau)x(10.2-\tau)\mathrm d\tau
= \int_{0}^{12} h(\tau)x(10.2-\tau)\mathrm d\tau
= \int_{7.2}^{12} h(\tau)x(10.2-\tau)\mathrm d\tau$$
(can you see why?)

So, there is not a simple answer to your question about the
limits of integration; you have to work them out for each
value of $t$ for which you want to evaluate the convolution
integral. But you don't have to do an infinite amount of
work!  There are Look Ma moments here and there.  Can you
combine two of the results above to write that for each $t$,
$3 \leq t \leq 9$, we have the useful pattern
$$y(t) = \int_{t-3}^{t+3} h(\tau)x(t-\tau)\mathrm d\tau$$
shown above?

I hope this will help you gain the understanding that
you need in order to complete your homework.


Get spectrogram matrix without specgram in MatLab

Yes, the easiest way to do this is using reshape and FFT, since reshape will give you a matrix and the FFT will operate along a dimension of the matrix. I've applied a 1024 rectangular window on a data set of 65536 points, sampling frequency Fs (1024), and signal frequency (100).

x = randn(1, 65536) + cos(2*pi*100/1024*(0:65536-1));
z = reshape(x, 1024, []);
y = 20*log10(abs(fft(z)));
imagesc(y)

This example doesn't have any overlap, but all you would need to do is shift z then concatenate it with z prior to taking the FFT. This all depends on your overlap percentage.

z = [z; circhift(z, [0 -1])];

This is not complete given the first column is now at the end, so you would have to omit that column, but you get the idea.


Why the difference come of 2D Fourier Transforms from two similar image?

Here is some code that will let you interactively make a mask and see the difference.

For the first image the dark edges and the bright centres are about the same width, so there is only one ring. For the second image, the dark edges are half the size of the bright centres, hence two rings.

I1 = imread('sefft1.jpg');
I2 = imread('sefft2.jpg');

I1g = rgb2gray(I1);
I2g = rgb2gray(I2);

I1gFFT = fft2(I1g);
I2gFFT = fft2(I2g);

imagesc(log(abs(fftshift(I1gFFT))))
imagesc(log(abs(fftshift(I2gFFT))))

% FFT coordinates
[rows,cols] = size(I1gFFT);

[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

% Centred versions
rc = fftshift(r);
thc = fftshift(th);

% Interactive selection of filter
imagesc(log(abs(fftshift(I1gFFT))));
title('Image 1: Press enter to select max and min freq');
pause;
title('select max and min freq');
[x,y] = ginput(2);
y = round(y)
x = round(x)
r1 = rc(y(1),x(1));
r2 = rc(y(2),x(2));
rmin = min(r1,r2);
rmax = max(r1,r2);
mask = (r > rmin) & (r < rmax);
maskc = (rc > rmin) & (rc < rmax);

imagesc(log(abs(fftshift(I1gFFT).*maskc))); title('Masked area'); pause;
imagesc(ifft2(I1gFFT.*mask)); title('Result of ifft of masked area'); pause;
imagesc(ifft2(I1gFFT.*(1-mask))); title('Result of ifft of non-masked area'); pause;

% Interactive selection of filter
imagesc(log(abs(fftshift(I2gFFT))));
title('Image 2: Press enter to select max and min freq');
pause;
title('select max and min freq')
[x,y] = ginput(2);
y = round(y)
x = round(x)
r1 = rc(y(1),x(1));
r2 = rc(y(2),x(2));
rmin = min(r1,r2);
rmax = max(r1,r2);
mask = (r > rmin) & (r < rmax);
maskc = (rc > rmin) & (rc < rmax);

imagesc(log(abs(fftshift(I2gFFT).*maskc))); title('Masked area'); pause;
imagesc(ifft2(I2gFFT.*mask)); title('Result of ifft of masked area'); pause;
imagesc(ifft2(I2gFFT.*(1-mask))); title('Result of ifft of non-masked area'); pause;


Image processing, recognizing a small feature in a larger image

What has been suggested to you is called template matching. If you don't mind re-using code (and working with C++), you can find an implementation in the OpenCV library.

A simpler alternative (that I would start with) would be to work with generalized Hough transform, but that will require some pre-processing first by detecting the edges in the image.


Your $C$ is here for the wavelet energy normalization: it is chosen such that:

$$ \int \psi(t)^2dt = 1\,,$$

the most  common framework for wavelets. It ensures that wavelet coefficients can more or less directly be interpreted in energy. The problem is that, as is, the wavelet is not admissible. If you take $\omega \simeq 0$, it reduces to a Gaussian, and does not  oscillate enough to have zero-mean.

However, when $\omega$ becomes "big enough", some  say $\omega>5$ (a common choice is $\omega = \pi\sqrt{2/\log2} $), one usually considers that it is "numerically admissible" (with high precision).

A more proper version (with $\omega_0$ instead of $\omega$ to avoid strange interpretations) is:

$$ \psi_{\omega_0} = \frac{C_{\omega_0}}{\pi^{1/4}} e^{-t^2/2}(e^{i\omega_0 t} -e^{-\omega_0^2/2})$$

and then you get:

$$ 1/C_{\omega_0}^2 = 1+e^{-\omega_0^2} -2e^{-3 \omega_0^2/4} \,,$$

which is very close to $1$ when $\omega_0$ is big enough.


Getting frequency content at different times from discrete wavelet transform coeffs

Question 1: I think you want to investigate centfreq &  scal2freq (Scale to frequency):

  FREQ = centfrq('wname') returns the center frequency in hertz of the
  wavelet function

  F = scal2frq(A,'wname',DELTA) returns the pseudo-frequencies
  corresponding to the scales given by A and the wavelet function
  'wname' (see wavefun for more information) and the sampling period
  DELTA.

The intuition that I use is to think about your question is: The wavelet transform is a hierarchical (tree-like) decomposition of an input signal. I would not say it is in the "time-domain" as is commonly understood. I think of the DWT output via the duality: Wavelets <---> Filterbanks.

Imagine a set of 8 bandpass filters tuned to filter different ranges from low to high on the frequency spectrum. You know which frequency ranges correspond to which filters, so you know how much signal lives in each of the 8 frequency ranges based on the bandpass filter output. That's a basic frequency analysis. You wouldn't need to compute an FFT at this point. And with the successive QMF filtering, then decimation of your signal at each level, your DWT coefficients may only be one-sample (or have a very coarse resolution), so an FFT wouldn't make any sense.

There's a good explanation @ https://math.stackexchange.com/questions/28581/which-time-frequency-coefficients-does-the-wavelet-transform-compute

Question 2: A scalogram (Scalogram (and related nomenclatures) for DWT? ) is good for visualizing your transformed signal, but I think it would be overkill just to compute which frequency range your wavelet coefficients correspond to.


Removal of noise from vibration recordings

If you want to detect the peak value of the time-domain signal, it is important not to change phase of each frequency signals. Otherwise, the peak value will be affected by it and will slightly change.

I think you already have recorded signal and you can process it after the recording.
If so, it is better to use zero-phased filters, like this function in matlab.

http://uk.mathworks.com/help/signal/ref/filtfilt.html

In this case, the needed function is to reduce the spike, which contains so much frequencies.
In the help page I have shown before, there is an example of filtering this kind of spikes from the relatively low-frequency signal.
I think it's what you need.

Edit

I recommend Butterworth low pass filter for its flat gain characteristics. The gain ripple of the Chebtshev would affect badly. The phase change will be corrected by using the filter backward, which the filtfit does. As you imagine, the effect of the phase change would be very little in this case, but it cannot be guaranteed.
Just for your information, Bessel would be a candidate when you do the filtering online.

If you still doubt about the effect of the phase shift (not the time-shift) by the filter, you can imagine the time-streched pulse and impulse signal in the image and code below.
These signals are almost the same but only the phase are different.
Besides the amount of the effect, the phase shift does affect to the amplitude of the signal.

I recomend you to use zero-phase filter if you want to measure the actual amplitude precisely.

clear all;
close all;

% generate signal with the same phase
N = 2^15; % number of samples for this signal
for k = 1:N/2;
    S(k+1) = exp(0);
end
S(N/2+2: N) = conj(S(N/2: -1 :2));

signal1 = real(ifft(S));

% generate signal with different phase from above
% this signal is called time-streched pulse.
N = 2^15;
for k = 1:N/2;
    S(k+1) = exp(-1j*2*pi*N/2*(k/N)^2);
end
S(N/2+2: N) = conj(S(N/2: -1 :2));

signal2 = real(ifft(S));

% and let's see them.
subplot(211)
plot(signal1)
xlim([-1e4, length(signal1)])
xlabel('time')
ylabel('all the same phase')
subplot(212)
plot(signal2)
xlim([-1e4, length(signal1)])
xlabel('time')
ylabel('only the phases are different form above')


Is there a public repository of labeled sound files coming from an industrial area?

Freesound is a repository of sound files categorised by user-defined tags. Through those, you can spot what other users have labeled industrial sounds but be prepared that some of those might be artificially created (so, not really a recording of an engine but something that sounds like an engine).

In the same category is also Soundsnap.

The other popular choice would of course be Soundcloud but there, you will have to do a bit of digging to find sounds that suit your use case because the majority of sound clips there are about music (so, "industrial rock" for example).

In all of these, you can curate a subset of samples that seem to have a strong relationship between their content and the tags, to then use the tag as your class variable if you are looking to download a large dataset.

Hope this helps.


How Does Mean Centering Affect the Result of Using SVD to Compress Images?

I will try to illustrate why it is important to remove the Mean from data when doing the PCA (SVD is the tool to basically do the PCA approach to dimensionality reduction -> Compression).

If the data isn't centered you will have one (Most probably the most dominant) component (Direction) which will point to the center of the data.
Since all directions must be orthogonal, it means the next directions, which should be correlated to the data itself will have harsh constraint.

You want the directions to be chosen according to the data layout and not according to its position in space.


GCC-PHAT always peaks at zero when estimating delay offset

This code seems to work (included also below).
Pythonizing the first example link you include seems to do the right thing with that code.
Fs = 8000
dt = 1/Fs #0.125e-3
f1 = 100
tdelay = 0.625e-03 # try different values
t3 = np.linspace(0.0,1.0,Fs)
x3 = np.cos(2*np.pi*f1*t3);
x4 = np.cos(2*np.pi*f1*(t3-tdelay));
gcc_phat(x4,x3)

yields:

(0.375,  array([-0.02514735, -0.02196879, -0.01847165, ..., -0.0303065
,
-0.02794452, -0.02514735]))

(which gives one minus the actual delay).
After playing with your wav file (thanks!), I think the problem is that your recording equipment has too much correlated noise across all four channels.  That's giving a GCC-PHAT peak at zero.
If I try to drown out the noise by adding even more (uncorrelated) noise, then I get a better result: -5 samples vs -4 samples from XCORR.

The four plots are:

The section of one channel that I'm looking at.

The GCC-PHAT central part when uncorrelated Gaussian noise is added to the data.

The GCC-PHAT central part when no noise is added to the data.

The XCORR result of the noiseless data.

Code to do this is below and on GitHub.  I originally though it might be possible to just filter the data, but that doesn't help at all.
import numpy as np
import statistics
from scipy.io import wavfile
import matplotlib.pyplot as plt
from IPython.display import Audio
from scipy.signal import kaiserord, lfilter, firwin, freqz

samplerate, data = wavfile.read("Q69905.wav",'rb')
samples = np.arange(20000,30000)

fir_filter = firwin(123,0.75)
channel_1 = lfilter(fir_filter,1, data[samples,0])
channel_2 = lfilter(fir_filter,1, data[samples,1])
channel_3 = lfilter(fir_filter,1, data[samples,2])
channel_4 = lfilter(fir_filter,1, data[samples,3])

noise_1 = np.random.normal(0,1000,len(channel_1))
noise_2 = np.random.normal(0,1000,len(channel_2))
noise_3 = np.random.normal(0,1000,len(channel_3))
noise_4 = np.random.normal(0,1000,len(channel_4))

print([statistics.mean(data[:,0]),  statistics.mean(data[:,1]), statistics.mean(data[:,2]), statistics.mean(data[:,3])])

# delay, gcc = gcc_phat(data[samples,0].astype(float)+10, data[samples,2].astype(float)+12, interp=1)
delay, gcc = gcc_phat(channel_1 + noise_1, channel_3 + noise_3, interp=1)
delay_no_noise, gcc_no_noise = gcc_phat(channel_1 , channel_3 , interp=1)

plt.figure(figsize=(20,30))
plt.subplot(4, 1, 1)
plt.plot(data[samples,0])
plt.subplot(4, 1, 2)
plt.plot(np.arange(-10,10),gcc[9990:10010],'.') # [9950:10050]
plt.subplot(4, 1, 3)
plt.plot(np.arange(-10,10),gcc_no_noise[9990:10010],'.') # [9950:10050]
plt.subplot(4, 1, 4)
lags,c, line, b = plt.xcorr(channel_1,channel_3)
plt.plot(lags,c,color='r')

print('GCC-PHAT: ' + str(delay))
print('XCORR: ' +  str(lags[np.argmax(c)]))
Audio(channel_1 + noise_1, rate=44100)

Code From GitHub
"""
 Estimate time delay using GCC-PHAT
 Copyright (c) 2017 Yihui Xiong
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
     http://www.apache.org/licenses/LICENSE-2.0
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
"""

import numpy as np

def gcc_phat(sig, refsig, fs=1, max_tau=None, interp=16):
    '''
    This function computes the offset between the signal sig and the reference signal refsig
    using the Generalized Cross Correlation - Phase Transform (GCC-PHAT)method.
    '''

    # make sure the length for the FFT is larger or equal than len(sig) + len(refsig)
    n = sig.shape[0] + refsig.shape[0]

    # Generalized Cross Correlation Phase Transform
    SIG = np.fft.rfft(sig, n=n)
    REFSIG = np.fft.rfft(refsig, n=n)
    R = SIG * np.conj(REFSIG)

    cc = np.fft.irfft(R / np.abs(R), n=(interp * n))

    max_shift = int(interp * n / 2)
    if max_tau:
        max_shift = np.minimum(int(interp * fs * max_tau), max_shift)

    cc = np.concatenate((cc[-max_shift:], cc[:max_shift+1]))

    # find max cross correlation index
    shift = np.argmax(np.abs(cc)) - max_shift

    tau = shift / float(interp * fs)

    return tau, cc

def main():

    refsig = np.linspace(1, 10, 10)

    for i in range(0, 10):
        sig = np.concatenate((np.linspace(0, 0, i), refsig, np.linspace(0, 0, 10 - i)))
        offset, _ = gcc_phat(sig, refsig)
        print(offset)

if __name__ == "__main__":
    main()


Image Denoising with Better Edge Preservation

Sounds like you want to denoise and preserve edges. Have you considered nonlocal means? There's some GPL'd C++ code along with a brief writeup of the algorithm by the original authors here: http://www.ipol.im/pub/algo/bcm_non_local_means_denoising/

One caveat, nonlocal means is very slow and the output can be sensitive to the implementation you have. You may also consider ROF minimzation as it's fast and does a good job of preserving edges. Here's some matlab code that does it: http://www.stanford.edu/~tagoldst/Tom_Goldstein/Split_Bregman.html


Ok I give the answer by myself, but it’s only half of answer. As I know we can rotate $ 90^o $ the phase by multiplying the wave by $ j $ - imaginary number. And to rotate it to $ 180^o $ we need multiply by $ j^2 $ (or just $ -1 $ which is the same). To rotate $ 270^o $ we need to multiply by $ j^3 $, and so on...

But if we want rotate something between? For example $ 45 ^o $ or $ 30^o $. What we need to multiply by?


What filter is better with -1/3 2/3 -1/3 or with -1/4 1/2 -1/4 coefficients?

Both filters can be re-written in the form:
\begin{equation}
H(z) = \alpha ( 1 -2z^{-1} +z^{-2} )
\end{equation}
Clearly, in the first case, $\alpha=-\frac{1}{3}$ and in the second $\alpha=-\frac{1}{4}$.

In other words, the shape of the filter's frequency response is the same, the only difference is the scaling.

I can think of two reasons of why the second filter may be preferred.

In the case of the second filter, the maximum gain of the filter (occuring at the Nyquist frequency) is normalized to $1$, whereas the first filter has a maximum gain of $\frac{4}{3}$.
The other more compelling reason however, is probably that the coefficients of the second filter can be expressed as powers of two ($\frac{1}{4}=2^{-2}$, $\frac{1}{2}=2^{-1}$) and therefore multiplication by these coefficients can be quickly implemented as bitwise operations on a computer or digital signal processor.


The Effect of a Kernel on an Image By Filtering

Convolution with following kernel returns the signal unchanged.

 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     1     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0

If the '1' is shifted to the right one position, the convolution output will be shifted to the right one position as well.

If the '1' is a '2' the convolution output will be double (scaled by 2).

If there is more than one non-zero number in the kernel, the convolution output will be the addition of each scaled and shifted copy.

For your kernel

 0     0     0     0     1     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     0
 0     0     0     0     0     0     1
 0     0     0     0     0     0     0
 1     0     0     0     0     0     0

the output is the addition of three copies of the signal, shifted by (1,3), (3,-1) and (-3,-3) and scaled by 1. I wouldn't call it blurring though.


smooth noisy irregularly spaced data containing peaks

The idea of cascading a non-linear then a linear filter is sound in the presence of spikes.

Here are some observations which, if correct, could be exploited for better processing:

Data is circularly arranged, and in fact 2D
Spikes point outward,
The remaining noise seems bounded (the thick blue ring),
The angular sampling is dense, and could be affected by noise.

Some suggestions:

Choose a novel angular sampling $\delta \phi$ that could be used for several datasets. It could be parametrized by a modelling of histograms with a one-sided exponential law, like $a\exp(-\lambda t)$, $t\ge 0$.
Choose an angular width $\Delta \phi$, that will be used as the span of a 2D filter.

Their combination would define a series of overlapping windows on which you would process 2D point clouds.

Inside each angular region defined by $[n\delta \phi - \Delta \phi/2, n\delta \phi + \Delta \phi/2]$, you can for instance:

define a nonlinear filter to eliminate strong spikes: for instance a 2D median, or a 2D weighted median (A fast implementation of two-dimensional weighted median filters). If you use the outward direction of spikes, an order filter, or a rejection of the top $x$% of the radius could be used. The idea behind weighting is to penalize the points that are the farer away from the angular sampling. A bilateral filter for point clouds could be an option too.

Introducing a weighting in a median filter can bridge the gap between pure median and smoothing filters, and kill two birds with one stone: remove spike, reduce noise. As a result, this could suffice. If not, you can then:

perform a 2D linear filtering on the remaining points, in the shape of a center of mass allocation.

One advantage with resampling to a regular angular grid is potentially to allow easiest further processing, comparisons, etc. And using 2D filters, you naturally follow the circular symmetry of your data.


How to get the data (pointwise frequency, amplitude and phase) in a FFT performed on ImageJ?

Here is a very practical illustration. The credit goes to @Herbie and @emartini.

NOTE: The complex window has a slide bar at the bottom, which has to be drawn to the left. Also, it is necessary to click twice on the FFT window without moving the cursor to get the three output values.

The output is concordant with Matlab with adjustments nicely described by @SleuthEye here and here.

Here is an example:

with ImageJ:

Frequency = 10.24 pixels/cycle (25 cycles)
Theta (direction) = 16.26 degrees
Real part = -1.255
Imaginary part = 10.142
Phase = arctan(10.142 / -1.255) = -82.95 degrees
Magnitude = sqrt(10.142^2 + 1.255^2) = 10.2194

and with Matlab:

Im = imread('circ.png');
pkg load image
Im = rgb2gray(Im);

A = fft2(double(Im) / 255);
Ashifted = fftshift(A);
Ashifted(121,153)

i = 121;
j = 153;

center = size(A) / 2 + 1;
dx = (j - center(2)) / size(A,2);
dy = (center(1) - i - 1) / size(A,1);

direction = (atan2(dy, dx))
dir_degrees = direction * (360 / (2*pi))
frequency = 1 /sqrt(dx*dx + dy*dy)

The output:

ans =  -1.2553 + 10.1425i
direction =  0.28379
dir_degrees =  16.260
frequency =  10.240


You can use the Normalized Cross Correlation for that.
http://en.wikipedia.org/wiki/Correlation_and_dependence

Basically, representing each recorded sound as a vector, this gives the angle between them.

Another approach is dealing with features of the sound instead of the signal itself.
You can start by reading how Shazam works:

https://www.google.com/webhp?sourceid=chrome-instant&ion=1&espv=2&ie=UTF-8#q=how%20shazam%20works

http://www.reddit.com/r/askscience/comments/1mdb08/how_do_apps_like_soundhound_and_shazam_work_is/

http://gizmodo.com/5647458/how-shazam-works-to-identify-nearly-every-song-you-throw-at-it

http://www.slate.com/articles/technology/technology/2009/10/that_tune_named.2.html

http://en.wikipedia.org/wiki/SoundHound

http://www.quora.com/How-does-SoundHound-work

http://everythingelsematterstoo.blogspot.co.il/2010/11/how-shazam-works.html


How to implement a high pass filter digitally to remove the dc offset from sensor data

The filter you desire is often called a “DC removal” or a “DC cancelation” filter. Unfortunately there’s no way to design a useful IIR DC removal filter that has no “settling time”. Your designed filter’s “transition region” is the lower passband frequency (50 Hz) minus the upper stopband frequency (10 Hz). So your designed filter’s transition region is 40 Hz.
If you design a “sharper” DC removal filter having a transition region of 20 Hz then that filter will have a longer settling time. If you design a DC removal filter having a wider transition region of 100 Hz (a lower “performance” filter) then that filter will have a shorter settling time than the filter you designed. So tryingengineer, you’re faced with a trade-off; the higher the performance of a filter (narrower transition region) the longer will be its settling time.
By the way, IIR filters have nonlinear phase responses. If you seek a computationally-efficient  DC removal filter that has a linear phase response, see the material at:
https://www.dsprelated.com/showarticle/58.php


Matrix form of 2D-DFT for a vectorized image

In matrix form, we can write the 2-D DFT of an Image $\mathbf X$ via $\hat{\mathbf{X}} = \mathbf F \mathbf X \mathbf F^{\rm T}$.

This expression can be vectorized with the help of the Kronecker product. Using [*], we obtain $$\hat{{\mathbf{x}}} = {\rm vec}\{\hat{\mathbf{X}}\} = (\mathbf F \otimes \mathbf F) \cdot \mathbf{x}.$$

This shows that the matrix you need is $\mathbf F \otimes \mathbf F$.


Sequential Non Linear Least Squares Problem

I have the the following non-linear function [...]

This would be a bimodal distribution with two different means and deviations ($\mu_1, \sigma_1, \mu_2, \sigma_2$), two sub population sizes $A_1, A_2$ and an extra parameter $t$ that modulates $A_1, A_2$.

  Here the data points DO NOT ARRIVE randomly from the model. They arrive in the order as given by the model. Namely, the incoming points are first around zero level, then they increase, make a peak, and then decrease, and the same happens for the second exponent.

What you are describing is plain simple curve fitting of some curve to data points as opposed to fitting data to a distribution.

The typical fitting scenario involves the use of Newton's method to find those parameters that minimise the value of the error function and here it is possible to differentiate the sum of the Gaussians with respect to the parameters (which is a common requirement for both linear and non-linear methods).

Here are a few practical examples: One, two, three (different platform).

If you can use details of your problem to fix the starting points of the fit, it will help immensely with convergence. There are some cases for example where the method converges but when you look at the fit, one curve has simply absorbed the other.

Can it be done with an online algorithm? No. Not in this case.

In general, there are online algorithms that will "learn" the coefficients of a model in a "sample-by-sample" basis. But this refers to working out parameters of a model that minimise an error function. Evaluating the error, implies a full "snapshot" of what is being fitted (i.e. all of the data to obtain one estimate of $s$).

Here, it is impossible to predict what sample comes next and from which possible realisation of the 6 parameter model. But there might be a thin region of cases where you could make an inference provided that you can apply additional constraints.

You might see this problem as a combination of a tracking and a classification problem where you can say "I know how the curve is supposed to evolve and all I have to do is track its evolution and raise a flag when the uncertainty falls below a level", but then again you don't really know which 6 parameter realisation to track. This is what you are trying to infer.

It is impossible to "guess the future" unless there are certain constraints.

What do the first few "... incoming points are first around zero level,..." tell you about the value of $A_1,\mu_1, \sigma_1, A_2, \mu_2, \sigma_2$? Nothing. You can still attempt a "fit" (batch or online) but you will get back nonsensical results and a huge error. So, we have to wait for more data. But then again, if we wait long enough to be sure we have one of the curves "in", then, we probably can make a quick guess about the $A$s provided that we know what $A_1+A_2$ is in general (another constraint). But still, this leaves us with no way of guessing what the $\mu_2, \sigma_2$ would be. We still have to wait for more data.

If the problem you are dealing with can already inform you about nonsensical values for the 6 parameters or at least the range they are supposed to be in or some sort of relationship between the two triplets of parameters, or anything else that could be used to reduce your search space, then, you might be able to pre-calculate a very large number of realisations of $s$ for all the different possible 6 parameter combinations and then track the evolution of each one as a new sample comes in. Which is basically, nearest neighbour classification.

But still, you would probably need to wait almost for the whole curve to come "in" to be able to make a certain decision.

If you start tracking a single Gaussian curve with your eye, you cannot really tell what $A, \mu, \sigma$ are until the trace starts curving to begin forming the peak. Add to this the noise and you have more than one possible curves that could be passing from the same data. So, you probably have to wait for at least $\frac{3}{4}$ of the data to come "in", before you can make a valid estimate.

Hope this helps.


What is the type of each 4 different noise signals that are shown below?

Let me start with Noise 4. Noise 4 is "Zero-mean White Noise".  White noise means the signal power is distributed across all frequencies evenly, which is what you see in the FFT result. From the histogram of the time signal we see that the signal is almost uniformly distributed, so close to a uniform white noise process. Although it appears closer to the sum of two independent Gaussian distributions, one with a positive non-zero mean and the other with a negative non-zero-mean, approximately +/- 0.25. (The result would be a zero-mean white noise random process).

Noise 1 is non-zero mean white noise, as evidenced by the spike at $F=0$ with a similar white noise spectrum everywhere else.  From the histogram we see that it appears to be the sum of a zero-mean Gaussian White Noise (AWGN) process with a non-zero mean random process. The power level in the non-zero mean process is so small compared to the zero-mean AWGN process that we are unable to really discern further properties of that non-zero part that is added except for the mean which appears as a strong component in the FFT. Also from the histogram it looks like there really is very little signal component between the two, so the time domain plot appears as it does simply because you are connecting the dots as a line plot rather than plotting just the points (try using plot(time, noise, '.') to see just the waveform points. I suspect that it is a Gaussian distribution with a mean of 10 that is sparse compared to the non-mean noise waveform.

Noise 2 has been low pass filtered, and looks like it may have possibly been a white noise process that has been filtered by adding to itself after having been delayed by one sample. This could be called Pink Noise since the low frequency noise dominates (drawing an analogy from the visible light spectrum). From the histogram we see that this is well approximated as a Gaussian random process. So in this case it would be band-limited Gaussian noise.

Finally Noise 3 has a negative DC offset. Because of the scale, we cannot see if the noise process is pink or white (would need to zoom in on the FFT to observe that.) However from the histogram we see that it is a Gaussian distribution with non-zero mean.


Null autocorrelation function and stationary

If E[X(t)] is constant and RX(t+τ,t)=0 can I say the process is WSS?
Can I say RX(t+τ,t)=0=RX(τ) and therefore is WSS?

Two times the same question.
It fulfills the definition (as you noticed yourself), so why even ask? yes.

What does RX(t+τ,t)=0 mean?

It means the process is uncorrelated ($R_X(\tau \ne 0)=0$), but also that its variance is 0 ($R_X(0) = \sigma_X^2 =0$). There's only one process that fulfills that:
$X(t) = 0$
not a really exciting process, is it?


Formally, the quantization is a mapping from a partitioned interval $[-A, A] = \cup_{k=1}^N I_k$ to the set of integers $\{1, \ldots, k\}$. The dequantization is the mapping from the quantization indices to a set of values $x_k$ chosen from each $I_k$.

When the quantizer used is uniform without any scaling (that is to say, when it is a simple rounding), the dequantizer is the identity. But when the quantizer is non-uniform, the integer passed into the entropy coder is not just a "rounded" value of the coefficient - an exponent, scaling factor, or any other monotonous transformation might have been applied to it. The dequantizer will have to reverse all the reversible steps in the quantification operation.

For example, in the AAC codec, the quantized value sent to the entropy coder is:

$$Q(x_i) = \text{sgn}(x) \left[\left(\frac{|x_i|}{s_i}\right)^\frac{3}{4} + 0.4054\right]$$

So the dequantizer will have to apply a $\frac{4}{3}$ exponent and multiply by the scale factor $s_i$ to recover the value.


Confusion regarding plot of linear convolution vs fast convolution via FFT

The variable NI in the book is just a typo, it should be N instead (not N*L as in your code). Apart from that, remember that book was written about $25$ years ago, and the code was run on a $33$ MHz $486$ PC. So in order to see some effects on today's computers, you should crank up the value of L.
I've modified the code a bit (see below). Now the figure illustrates the expected result.

L = 200:500:10000;
K = length(L);
conv_time = zeros(1,K); fft_time = zeros(1,K);
Nav = 10;

for k = 1:K
    tc = 0; tf=0;
    Lk = L(k);
    N = 2*Lk-1; nu = ceil(log10(N)/log10(2)); N = 2^nu;
    for I=1:Nav
        h = randn(1,Lk); x = rand(1,Lk);
        t0 = clock; y1 = conv(h,x); t1=etime(clock,t0); tc = tc+t1;
        t0 = clock; y2 = ifft(fft(h,N).*fft(x,N)); t2=etime(clock,t0);
        tf = tf+t2;
    end
    conv_time(k)=tc/Nav; fft_time(k)=tf/Nav;
end

plot(L,conv_time,L,fft_time)


Is step detection the correct approach to this problem? what if not?

First off, I would start by applying a smoothing/averaging filter to your position data. This will get rid of a lot of the noise from your stationary motion.

To answer your first question, step detection is the way to go, but it depends on how you're detecting the steps. The best way in my mind would be to look at the derivative, the speed, and detect periods of non-zero speed. If you've smoothed out your position data, this should result in periods of time where the animal is moving in one direction, rather than back and forth.

This approach makes two dimensional analysis fairly simple. Once you've calculated your x and y speeds, you can combine them to get an absolute speed, independent of direction. You can use the same detection approach to find the periods of motion, which should give you the answer you're looking for.


Carrier Frequency Offset Correction Limitations

There is no fundamental limit but based on what frequency discrimination approach you are using but you will be trading sensitivity for range. To improve sensitivity you can use a cascade of discriminators to "tune in" on your signal of interest.  A trivial example of this is a case with a carrier at 2.4 GHz; regardless of what carrier frequency we start with as our initial "guess", we could band-bass filter the signal of interest with a fixed filter at 2.4 GHz (to eliminate content at other frequencies that would disrupt the discriminator) and then use a succession of frequency discriminators (which simply provide an output voltage proportional to the frequency of the input) to provide higher and higher accuracies of what the carrier is, which can ultimately be within the capture range of whatever you are using for OFDM. For example, this link provides a typical discriminator for OFDM:

Carrier frequency offset (CFO) estimation using cyclic prefix

And this link provides further details on carrier discrimination techniques in general:

What exactly is a 90 degree phase shift of a digital signal in FM demodulation appraoches?


What is the autocorrelation of a Dirac pulse?

Well, by definition of the $\delta$ distribution, you have:

$\int_{-\infty}^{\infty} f(t) \delta(t-T)\, \textrm{d}t = f(T)$

The autocorrelation of a function $g(t)$ can be computed via:

$\int_{-\infty}^{\infty} g^{*}(t)g(t + \tau)\, \textrm{d}t$, with $g^*$ as the complex conjugate of $g$. Since $\delta(t)$ is real-valued, this is conjugation can be skipped. So you are left with:

$\int_{-\infty}^{\infty} \delta^{*}(t)\delta(t + \tau)\, \textrm{d}t = \int_{-\infty}^{\infty} \delta(t)\delta(t + \tau)\, \textrm{d}t = \delta(-\tau)$.

The first = sign comes from the autocorrelation of the real-valued $\delta$, the second from the definition of the $\delta$-distribution.

So, the autocorrelation function of the $\delta$-distribution is the distribution itself. A eigenfunction of the autocorrelation function, so to say ;)

Think about it, this does make sense: the only perfect match is achieved with no time shift, ie at $\tau = 0$. All other shifts would end up with one of the arguments of the $\delta$ being different from 0, hence with the $\delta$-function being 0 there.

BTW: $\delta(-\tau) = \delta(\tau)$, since the function/distribution is symmetric.


There are a few different ways to proceed. They involve doing some algebra but it is more an issue of tedium. The notation of your analog filters is a little non standard, essentially your $p$ is $s$ in most books.

Traditionally, phase wasn't considered important in hearing so standards tended to be specified by magnitude, so there is some ambiguity with respect to the digital implementation's phase response. You have 2 basic simple choices, a digital IIR filter or FIR filter. FIR is likely to be easier but easy and better are typically not the same. FIR can be direct or FFT based.  IIR has a number of choices such as biquads or direct. Floating point math has its own easy/better issues.

One path way is to focus on the terms inside the magnitude brackets, i.e. the complex analog transfer functions. The other path is to focus on the analog magnitude.

So, given the analog transfer functions, you can either use a table of Laplace transforms and deduce the impulse response of each component you listed, or use the bilinear transform to express each component that is in $s$ in terms of $z$.

If you derive the impulse responses, you can sample points of the impulse response and provided a generally exponential decay, use those samples as coefficients of a FIR filter of sufficient length. The FIR filters are then applied sequentially. One nice property of ideal LTI filters, the order that you apply the filters are ideally interchangeable.

If you choose the bilinear transform and do the algebra to derive the discrete time numerators and denominators, you need to be aware that the bilinear transform squeezes the response as it approaches the Nyquist frequency. One needs to plot the discrete time response to asses how the corresponding analog response is distorted. If some modest tweaking results in an acceptable response then, the numerators and denominators need to be factored so that they are compatible with the IIR implementation, such as a cascade of biquads.

If you go with the magnitudes, the task becomes a direct sampling from the analog frequency domain to the discrete time periodic frequency domain. One can choose a linear phase response, and appropriately impose symmetry and define the frequency domain FIR filter. This might not sound like an analog implementation.

The essential challenge is mapping a frequency response in $-\infty < \Omega < \infty$ to $-\pi \le \omega \le \pi$

There are other approaches like state variables, so this is not the entire story.


Finding circle in binary image (R, not matlab)

On page 4 of this document describing the R PET package, there's documentation of the Hough transformation. That's probably what you want.


UPDATE: On choosing the size of $R_{xx}$: If you suspect there are $M$ sources, then you need your autocorrelation matrix ($R_{xx}$) to be made up of M+1 unique (i.e. non-overlapping) snapshots.  For this reason, the snapshots technically have to be less than len(x) / M+1.  Estimating $R_{xx}$ is not trivial, and there are some pitfalls for certain signals and certain selections of N (size of autocorr matrix) if the segments or slices are non-overlapping.  This is why I typically use as many segments of x as possible when estimating $R_{xx}$

First, you need an autocorrelation matrix.  One of the more fool-proof ways to do that, is to add up successive outer products of "slices" of the input vector where a slice is x[idx:idx+N].

# x contains input vector , usually a time-series
N = 20 # length of one side of the square autocorr matrix
Rxx = NxN matrix of zeros
foreach unique slice of x of length N:
    s = unique slice
    Rxx += outer product of (s,s) #outer(s,s) in Numpy
Rxx /= nsnap # divide by number of snapshots that Rxx is composed of

Second, calculate the eigenweights and eigenvectors of $R_{xx}$.  Save the "minimum" eigenvector. The minimum valued eigenweight and minimum eigenvector should correspond to the subspace of the noise.

w,v = eig(Rxx)
vmin = v[:,abs(w).argmin()]

At this point, you have $v_{min}$ and $e^H$ is simply

eh = exp( 1j*omega*arange(N) )

Where $PHD = \frac{1}{ | e^H v_{min} |^2}$ or

phd[omegaIdx] = 1.0 / ( abs(eh * vmin) ** 2 )

Do that for all omegas of interest.  The non-pseudo-code is in Python/Numpy, hence a * b is an element-by-element multiplication.


Fill matlab "ellip" function through transfer function magnitude response

The quickest way would be

[b,a] = ellip(8,1.5,60,[.2 .3]);

This designs a bandpass from .2 to .3 with 1.5 dB ripple and 60 db stop band attenuation. You don't get to pick the edges of the stop band since you already pre-selected the order. You can do one or the other but not both at the same time. However the 8th order filter is more than enough to meet your stop band requirements. In fact the low band stops at 0.195 and the high stop band starts at 0.307.

Another wrinkle is that the ripple actually only goes down, i.e. passband is not between -1.5dB and +1.5 dB but only between -1.5dB and 0 dB. In order to work aound this (if you care) you can double the ripple to 3 dB and than simply add 1.5 dB gain. You would also have to increase the stop band attenuation to 61.5dB to compensate for the overall gain change. The funcion $ellipord$ allows to calculate the required order to achieve a certain stop band, so the correctest way would be

%% Ellip filter design example
[n,wp] = ellipord([.2 .3],[.14 .36],3,61.5);
[b,a] = ellip(n,3,61.5,wp);
% add 1.5 dB of gain
b = b*10.^(1.5/20);
[H,F]= freqz(b,a,1000,2);
plot(F,10*log10(H.*conj(H)));
set(gca,'ylim',[-70,3]);
grid('on');

So it turns out a 5th order elliptic filter can do the job just fine.


In the multi-path scenario we have occurrence energy from (typically) all directions, and a channel could be a Rayleigh fading channel; no LOS component or Rician fading channel; a strong LOS component is present.

The figure you referred to contains a dominant component with a delay spread of 130 ns and with responses from other paths occurring with different delay spreads. Also, the key to the figure mentions an RMS delay spread of 65.9 ns, which can be justified by noticing the higher power distribution between 0 to 100 ns delay spread.

Take a look at the picture below to better understand the CIR for a multipath scenario. In this case, the CIRs from multiple antennas will be combined to form the CIR of the antenna array.


Why do we need multiple layers in each octave and multiple octaves in SIFT?

In music theory, an octave is an interval in frequency, from a frequency $f$ to frequency $2f$. For example "an octave higher" means "twice the frequency". Expressed as wavelength inversely proportional to frequency, $\lambda \propto \frac{1}{f}$, an octave would be the interval from a $\lambda$ to $\frac{1}{2}\lambda$. In the SIFT paper's notation, $\sigma$ is the scale, which is analogous to wavelength as it is proportional to the width of the Gaussian kernel.
Downsampling is just to reduce the image size so that there is less computation to do. From the paper (emphasis mine):

Once a complete octave has been processed, we resample the Gaussian
image that has twice the initial value of $σ$ (it will be 2 images from
the top of the stack) by taking every second pixel in each row and
column. The accuracy of sampling relative to $σ$ is no different than
for the start of the previous octave, while computation is greatly
reduced.

What enables to do this? A Gaussian function has another Gaussian function as its Fourier transform, and that decays faster than can be counteracted by any differential operator such as the Laplace operator. A finite-amplitude low-pass signal filtered by a Gaussian or a Laplacian of Gaussian or its approximation of scale $\sigma$ has in an approximate sense a band limit that is proportional to $\frac{1}{\sigma}$. By increasing the constant of proportionality, the mean square power of the part of the signal that is above the band limit can be made arbitrarily small. This allows to resample the filtered signal with a sampling frequency proportional to $\frac{1}{\sigma}$ while keeping the aliasing error arbitrarily small. A convenient approach is to not sample every signal filtered at scale $\sigma$ by sampling frequency $\frac{\sigma_0}{\sigma}$ with $\sigma_0$ some reference scale that is sampled at sampling frequency $1$, but instead by sampling frequency $2^{\operatorname{ceiling}\left(\log_2(\frac{\sigma_0}{\sigma})\right)}$:
Figure 1. Sampling frequency (red) is the same for all frequencies within an octave. The drawback is that for some scales the sampling frequency is higher than the actual requirement (blue) for that scale. $\sigma$ is the scale and $\sigma_0$ is a reference scale.
This bunches together values of $\sigma$ into a number of octaves such that within each octave the sampling frequency is a power of 2, ensuring that the location of each spatial sample conveniently coincides with a spatial sample one octave higher. The drawback is that the sampling frequency for all scales in an an octave must be chosen based on the lowest scale included in the octave, to keep the aliasing noise within acceptable limits for that scale. The higher scales in the octave will typically have less aliasing.


Using MFCC to an ANN Speech Recognition System

You should process the speech in blocks of a few milliseconds at a time. Each block will be the same length. Then you can apply your ANN to those features.


IIR filter implementation in Direct Form 2

It seems to be wrong. The first sample of the impulse response is always equal the b0 coefficient (b(1) in Matlab, so in this case  that should be 0.099825504845559299 which doesn't seem to be the case.

The correct answer can be checked with yref = filter(bs,a,x); and it's substantially different from your answer. You have shown that the fixed point version matches the floating point version, however the floating point version is wrong to start with.

Finally, Direct Form II is a horrible choice for fixed point filters. The transfer function between input and state variables is given by the pole-only transfer function. Even in this harmless example this is about 27 dB of gain. It's almost impossible to prevent clipping in Direct Form II. Use Direct Form I or Transposed Form II instead


Question about the derive of filter roll-off

For any arbitrary rational complex function $H(\omega) = \frac{N(\omega)}{D(\omega)}$,

$$
|H(\omega)| = \frac{|N(\omega)|}{|D(\omega)|}
$$

This follows straightforwardly from the properties of the polar (i.e. magnitude/phase) representation of complex numbers. So, for your example:

$$
|H(\omega)| = \frac{|1|}{|1 + j \frac{\omega}{\omega_c}|}
$$

The magnitude of an arbitrary complex number is:

$$
|a + jb| = \sqrt{a^2+b^2}
$$

So for your example:

$$
\begin{align}
|H(\omega)| &= \frac{1}{\sqrt{1^2 + \left(\frac{\omega}{\omega_c}\right)^2}} \\
&= \frac{1}{\sqrt{1 + \left(\frac{\omega}{\omega_c}\right)^2}}
\end{align}
$$


How to keep phase response of IIR filter?

To my knowledge, you can't linearize it arbitrarily, but that being said, you can obtain approximately linear phase by processing the output through a cascade of all-pass filters.

An all-pass filter used this way is usually called a phase compensator or a phase equalizer. A single pole-zero all-pass IIR filter has the following structure
$$ H(z) = \frac{p - z^{-1}}{1 - p z^{-1}} $$
where $p$ controls the location of the pole/zero. You can optimize a cascade of filters of this type to equalize your phase response using least squares or some other filter design technique.

Note that this technique reduces the primary advantage of choosing an IIR filter to begin with. Usually, the point in choosing an IIR filter over an FIR filter is to reduce computational complexity, so performing the phase equalization by cascading all-pass sections might result in a situation where using a linear phase FIR filter in the first place would have been more economical.


Possible to implement simple IIR filter with greater than 6 db/octave slope?

You need a filter with a higher filter order, which means that your filter memory must be longer. The current output sample $y[n]$ is computed as

$$y[n]=a_1y[n-1]+a_2y[n-2]+\ldots + a_Ny[n-N]+b_0x[n]+b_1x[n-1]+\ldots + b_Nx[n-N]$$

where $x[n]$ is the current input sample, $N$ is the filter order, and $a_i$ and $b_i$ are the filter coefficients. You can design the filter (i.e., compute the coefficients) using some tool like Matlab or Octave.


Proof of Paley-Wiener criterion for causality

This is nice question which I have recently encountered while studying ideal filters. Suppose transfer function $H(p)$ which is complex rational and holomorphic (or meromorphic?) in $p = \sigma + j\omega$ and satisfy Cauchy-Riemann conditions.

We can rewrite transfer function as $H(p) = e^{-\Gamma(p)}=10^{-\frac{1}{20}\Gamma(p)}$ with respect to nepers and decibels. Since exponetital function is holomorphic so the $\Gamma$ is. If we restrict ourselves to imaginary axis $p = j\omega$ we have

$\Gamma(j\omega) = A(j\omega) + jB(j\omega)$

where $A(j\omega)$ is an attentuation charateristics and $B(j\omega)$ is a phase characteristics.
In general we have identity $\ln z = \ln |z| + j\arg z$ for any $z \in \mathbb{C}$. So we can write

$A(j\omega) = \ln|H(j\omega)|$ and $B(j\omega)=\arg H(j\omega)=\arctan\frac{\Im\lbrace H(j\omega)\rbrace}{\Re\lbrace H(j\omega)\rbrace} $.

Then consider pair of Hilbert transforms given by integrals (Cauchy principal value are taken in each of the integrals)

\begin{equation}
A(\omega) = \frac{1}{\pi}\int_{-\infty}^{\infty}\frac{B(j\xi)}{\omega - \xi}\,\mathrm{d}\xi \quad\quad B(\omega) = -\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{A(j\xi)}{\omega - \xi}\,\mathrm{d}\xi
\end{equation}

Evaluating these integrals for module characteristic of ideal low pass filter - say let $A(\omega)$ is constant at $a$ from $-\omega_1$ to $\omega_1$ (module characteristics are symetric) gives
\begin{equation}
B(\omega) = - \frac{a}{\pi}\ln\frac{|\omega + \omega_1|}{|\omega - \omega_1|}
\end{equation}
Thus $B(\omega)$ has points in infinity at $\pm\omega_1$.

See picture bellow:

Similarly for linear phase characteristics $B(\omega) = a\omega$.

This is saying to us that we can't reach filter transfer function $H(p)$ with both ideally rectangular module characteristics and ideally linear phase characteristics at the same time.

We see that $B(\omega)$ with points at $\pm \infty$ can't realise any phase characteristics $\arg H(j\omega)$ of a given feasible filter. That's the reason why
\begin{equation}
-\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{A(j\xi)}{\omega - \xi}\,\mathrm{d}\xi < \infty
\end{equation}
must converge and so the
\begin{equation}
-\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{|A(j\xi)|}{\omega^2 + \xi^2}\,\mathrm{d}\xi < \infty
\end{equation}
must.


Implementation of wideband beamformer for planar array

You need the polar response of each individual array element. That is the transfer function of each element as a function of azimuth and elevation (and potentially distance if you want to operate in the near field). For many arrays the polar function will be the same for all elements, but it doesn't hurt to check.
To calculate the polar response of the whole array you need to do the following

Create a spatial grid on which to calculate the polar response (azimuth and elevation) for the array
For each point on the grid calculate distance, azimuth and elevation of this point AS SEEN BY THE INDIVIDUAL ELEMENT. This is in essence a transformation between two different spherical coordinate systems.
Look up the single element transfer function for this direction in the single-element polar data. This may require interpolation between measured directions
Apply this transfer function to the input signal of this individual element
Make sure to properly account for the distances. That's typically done by a factor exp(-j*w*(r-r0)), where r is the distance between grid point and element and r0 and convenient normalization distance. A good choice is the distance at which the single element polar was measured.
Some results over all elements
Repeat for all points on the target polar grid


what is the difference between QPSK and 4QAM?

There are different types of QPSK: Offset (staggerd) QPSK, $\pi/4$-QPSK, etc.

In its most known form, QPSK is identical to $4$-QAM or $4$-PSK. That is, the same constellation can be referred to by different names. So there is no difference.

This constellation can be seen in the picture bellow (from Wikipedia)

Similarly, $2$-QAM is actually the same as BPSK, as is also pointed out here.


Classifying this signal in dimension and number of channels

I don't feel like it's a very well-posed question, as multi-channel signals are sometimes represented as vector-valued, which would make them multidimensional. However, I assume the answer that is being fished for is:

This is a one-dimensional, two-channel setup. The two channels are the child's height and weight. Each signal is one-dimensional, measured across time. The sample rate of each signal is $1\text{ month}^{-1}$. You could state them mathematically as:

$$
h[n] \equiv \text{ height of child at month } n
$$
$$
w[n] \equiv \text{ weight of child at month } n
$$


I suppose that conversions from signed 2's complement to Q15 and vice versa are necessary.

Yes & No.  Q15 (signed) uses 2's complement, the bits are the same so no conversion is necessary. The difference between a Q15 and Q0 aren't the bit them selves but your interpretation of them. If you multiply two 16-bit numbers you get a 31-bit long result (one less than 32 because of a redundant sign bit). If you want to convert this to Q15 you take the top 16 bits, if you want Q0 you take the bottom 15 bits plus the sign bit.
The problem is that fixed point audio processing is difficult and 16-bit audio processing is extremely difficult. You need to constantly manage head room and optimize signal to noise ratio for each processing step. Fixed point filters require careful management of section ordering, section gain, section topology, noise shaping, rounding strategies, limit-cycle prevention, etc.
If you can afford it, I strongly recommend to do this in 32-bit floating point, which is much, much easier. If you can't, consider 32-bit or at least 24-bit fixed point. If that's not possible either, at the very least use double precision (fixed point) accumulators and filter state variables. If that's also not an option, you have a very rough project ahead of you.


Discrete Fourier transform - Norms of complex input signals and their transformation

I may consider them as L^2 then following theorem holds.

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

https://de.wikipedia.org/wiki/Satz_von_Plancherel

https://ms.mcmaster.ca/craig/craig-pde-chap5.pdf

Also keen to learn if there is any other theorems.


How important is it to use power of 2 when using FFT?

Modern FFT libraries, such as FFTW and Apple's Accelerate framework can do non-power-of-2 FFTs very efficiently, as long as all the prime divisors of the composite length are fairly small (2,3,5,etc.)

A power of 2 makes it simpler (about 1 page of source code) if you have to code your own FFT for some reason, or are otherwise constrained as to max program length (or FPGA gates, etc.)

For phase measurement, it might be easier to do an fftshift (pre-rotate the data by N/2) to reference FFT phase to the center of the data window, where the evenness/oddness ratio, and thus the phase won't change or alternate with bin number (for phase that is the same at the center of that data window) even for signals that are non-periodic in the FFT length, as you vary the length.


Determine the optimum receiver and the corresponding $P_{eM}$ for an AWGN channel

As has been referenced in the comments already, you're describing a pulse-amplitude modulation (PAM) signal constellation. The problem, as you've framed it, seems to suggest the AWGN vector channel (where symbols are described using discrete values $s_1$, $s_2$, and so on), in contrast to the waveform channel, where symbols are expressed using waveforms that are transmitted over some period of continuous time.

While in the real world, the waveform channel is often used as a communication medium, analysis of performance of digital communication systems is typically related back to the vector channel for ease of analysis. This document has a discussion of the distinction between the two channels and how the conversion between the two is modeled mathematically.

With that said, you can frame the problem rather simply. Assume that at each symbol time, your receiver observes a value $x$, which is defined as follows:

$$
x = s_i + n,
$$

where the value of $i$ indicates the symbol that was transmitted. $s_i$ refers to the value that would be expected at the receiver in the case of no noise (i.e. $s_i$ refers to the amplitude that is transmitted for symbol $i$). I will state without proof that the optimal decision rule (for minimum probability of bit error) for the case of equiprobable input symbols and white Gaussian noise is to pick the symbol value whose constellation position is closest to the observed point. This result follows from Bayesian decision theory.

The exact expression for the bit error rate will depend upon the coding of bits onto the symbols in your constellation. Borrowing a figure from an above comment:

These are examples of Gray coded constellations, such that adjacent symbol positions only differ by one bit. These codings are used to minimize the probability of bit error for multi-bit-per-symbol modulations. You might ask: why do you need to know the symbol encoding in order to calculate the system's error performance?

The answer might not be obvious: the overall bit error probability is a single number that is averaged across the ensemble of all possible transmitted symbols and observed channel noise. However, it should be clear from the above figure that the bit error probability must be different based upon which symbol is actually transmitted.

  Example: Assume the $M=8$ case illustrated above, with a transmitted symbol of $000$.

  The noise observed at the receiver $n$ can take on any value $< 0$ without having any chance of causing a bit error! $n$ could be $-10000000a$ (which is likely not very probable), but $000$ would still be the closest constellation position and therefore the receiver would make the right decision. We clearly couldn't say this for another transmitted symbol like $110$, however, so the bit error probability must be dependent upon the symbol that was actually transmitted.

Assuming a particular coding of bits onto the constellation (e.g. the illustrated Gray code), you can construct the error probability conditioned on the actually-transmitted symbol. I will tackle the $M=4$ case from above because it is simpler than the 8-point constellation, but the same principle applies. While the spacing between points isn't shown on the graph (because I'm too lazy to draw one), the points from left to right should be located at $-1.5a$, $-0.5a$, $0.5a$, and $1.5a$, respectively. Remember that the optimum receiver chooses the constellation point that is closest to the observed value.

$$
P(\text{error | 00 transmitted}) = \frac{1}{2} \\
(P(\text{choose 01 | 00 transmitted}) + \\
2 P(\text{choose 11 | 00 transmitted}) + \\
P(\text{choose 10 | 00 transmitted}) )
$$

So we've very straightforwardly stated what the probability of a bit error is for this case. If $00$ was transmitted and $01$ is chosen by the receiver, then you have one bit error; if $11$ is chosen by the receiver, you have 2 bit errors; and so on. Since we know the decision rule employed by the receiver, we can state this more directly in terms of the observation. Remember that if the symbol $00$ was transmitted, then the observation should have the form:

$$
x = s_i + n = -1.5a + n
$$

$$
P(\text{error | 00 transmitted}) = \frac{1}{2} \\
(P(-a < x < 0\ |\ x = -1.5a + n) + \\
2P(0 < x < a\ |\ x = -1.5a + n) + \\
P(x > a\ |\ x = -1.5a + n) )
$$

We can rearrange the above to be solely in terms of the noise term $n$:

$$
P(\text{error | 00 transmitted}) = \frac{1}{2} \\
\left(P(0.5a < n < 1.5a) + \\
2P(1.5a < n < 2.5a) + \\
P(x > 2.5a)\right)
$$

On the AWGN channel, $n$ is a sample from a zero-mean white Gaussian random process. Assuming the noise process has variance $\sigma^2$, the above probabilities can be written in terms of the Gaussian distribution's tail probability function, known as the Q function:

$$
P(\text{error | 00 transmitted}) = \frac{1}{2} \left(\left(Q\left(\frac{0.5a}{\sigma}\right) - Q\left(\frac{1.5a}{\sigma}\right)\right) + 2\left(Q\left(\frac{1.5a}{\sigma}\right) - Q\left(\frac{2.5a}{\sigma}\right)\right) + Q\left(\frac{2.5a}{\sigma}\right)\right)
$$

In most cases, the error probability will be dominated by the closest constellation points (the noise is unlikely to be large enough to cause you to skip multiple points over in the constellation, so the corresponding probabilities will be small compared to those for the nearby points), but for strict correctness, you would need to include all of the terms. You can then repeat this analysis for each constellation point, calculating what the probability of making all of the incorrect choices are. Then, to arrive at a single bit error probability number, you use the law of total probability:

$$
P_e = \sum_{i=1}^{M} P(\text{error | }s_i\text{ transmitted}) P(s_i\text{ transmitted})
$$

Since all of the symbols are equally probable,

$$
P(s_i\text{ transmitted}) = \frac{1}{M}
$$

$$
P_e = \frac{1}{M} \sum_{i=1}^{M} P(\text{error | }s_i\text{ transmitted})
$$

And that is how you would calculate the bit error performance for an arbitrary PAM constellation. Note that it is inherently a function of the signal-to-noise ratio, implied via the proxy values of $a$ (which is related to the received amplitudes, which is related to the received "energy per bit" $E_b$) and the noise variance $\sigma^2$ (which is also often dubbed $N_0$ in digital communications analysis). You could manipulate the above for given values of $M$, $a$, and $\sigma^2$ to arrive at a final result that is a function of the typical metric $\frac{E_b}{N_0}$. That is left as an exercise to the reader.


OK, so let's have a look at whether

$$\sum_{m=0}^{N-1} \sum_{n=0}^{N-1} |u(m,n)|^2 = \sum_{k=0}^{N-1}\sum_{l=0}^{N-1} |v(k,l)|^2\tag{1}$$

is equivalent to

$$U^2=V^2\tag{2}$$

Let's have a look at a simple $2 \times 2$ example:
$$
U = \left [
\begin{array}{cc}
u_a & u_b\\
u_c & u_d
\end{array}
\right ] \\
V = \left [
\begin{array}{cc}
v_a & v_b\\
v_c & v_d
\end{array}
\right ]
$$
so that
$$
U^2 = \left [
\begin{array}{cc}
u_a & u_b\\
u_c & u_d
\end{array}
\right ]
\left [
\begin{array}{cc}
u_a & u_b\\
u_c & u_d
\end{array}
\right ] = \left [
\begin{array}{cc}
u_a^2 + u_b u_c & u_au_b + u_bu_d\\
u_cu_a + u_d u_c & u_c u_b + u_d^2
\end{array}
\right ]
$$
and
$$
V^2 = \left [
\begin{array}{cc}
v_a^2 + v_b v_c & v_av_b + v_bv_d\\
v_cv_a + v_d v_c & v_c v_b + v_d^2
\end{array}
\right ]
$$

Then (1) becomes
$$
|u_a|^2 + |u_b|^2 + |u_c|^2 + |u_d|^2 =
|v_a|^2 + |v_b|^2 + |v_c|^2 + |v_d|^2
$$

The two are completely different.  Why do you think they are the same?


Find cut-off frequency of a low-pass filter for a given output signal

set $\theta_c = \pi$ and set $x[n] = y[n]$ and you're done.  and you're not dividing by zero anywhere.


Kalman Filtering with Unknown State Transition Matrix

Unfortunately, in Classical Kalman filter applications you need to have a guess of your state transition model governing the states. Classical choices for such models for time-varying systems are the random walk model, or its generalization, a first-order Markov model. If do not have a clue about which state model you can assume, you can try to extract a trend component from your state time series (via EMD or a moving average filter, for example), and plot it to have an idea of how the first-order moment of your states vary with time, for example. If it decays as a power law with $n$ (i.e. $E[x(n)] \approx \alpha^n$) until a final value, it may be an indication that a first-order Markov model with pole equal to $\alpha$ can serve as a state transition model.

Otherwise, you can treat the transition model for your states as a black box and consider techniques such as the Unscented Kalman Filter, for example.


The equation looks reasonable to me: The scaling by $\frac{1}{2\pi}$ is to have the result in units of normalized frequency $f$ instead normalized angular frequency $\omega$. What may be confusing is using the index $k$ from the autocorreation instead of $n$ since it would be the time domain variable for $R[k]$ and $w[k]$, while $k$ is often associated with the frequency index in the DFT. The OP is seeking the frequency variable to be in the exponent, in this case suggesting $n$, but instead the unit used is actually $\omega$ representing the fundamental normalized angular frequency (which for a DFT over $N$ samples would be given as the expected $2\pi n/N$ as in the OP's 2nd question).

In this case, there are $2N-1$ total samples consistent with a linear correlation over $N$ samples apparently given as:

$$R[k] = \sum_{n=0}^{N-1}x[n]x^*[n-k], \space\space\space k \in [-N-1, \ldots N-1]\tag{1}\label{1}$$

The Discrete Fourier Transform of the windowed autocorrelation function could then be (using a frequency index $\ell$ and assuming there are also $2N-1$ frequency samples):

$$X(\ell) = \sum_{k=N-1}^{N-1}w[k]R[k]e^{-j2\pi k \ell /(2N-1)}, \space\space\space \ell \in [-N-1, \ldots N-1]\tag{2}\label{2}$$

This will grow by a factor of $2N-1$, which we could normalize if desired by dividing by $2N-1$. As in any DFT, this is just a scaling and its used would be dictated by defining the units of the output.

Finally the normalized radian frequency (the fundamental frequency when $\ell = 1$) in this case is $\omega = 2\pi \ell/(2N-1)$, if substitute this in $\ref{2}$ we get a function of $\omega$:

$$X(\omega) = \sum_{k=N-1}^{N-1}w[k]R[k]e^{-j\omega k}\tag{3}\label{3}$$

Or if instead if we want a function of $f$, given $\omega = 2\pi f$:

$$X(f) = \frac{1}{2\pi}\sum_{k=N-1}^{N-1}w[k]R[k]e^{-j\omega k}\tag{4}\label{4}$$

I was not able to locate an actual copy of Blackman and Tukey's 1958 paper (which predated the FFT) but did find this additional variation from "Noninvasive Instrumentation in Medical Diagnosis" by Robert B. Northrop. In this variant, I would think that either the exponent should be divided by $(N+1)$ or the summation and number of output frequencies should be done over $N$ samples:

The OP asked in comments if this is best practice for estimating power spectral density. The power spectral density is given as the Fourier Transform of the auto-correlation function, and doing a linear correlation (which can be done by zero-padding the time domain sequence and doing a circular correlation) and windowing the resulting time domain sequence prior to computing a DFT are indeed best practice for estimating the Fourier Transform. The linear correlation avoids time domain aliasing, and the windowing reduces spectral leakage error. So this is standard practice.

One precaution not yet mentioned with windowing and spectral estimation: the equivalent noise bandwidth of a rectangular window is 1 DFT bin, which means when the spectrum in generally equally distributed the total power is equivalent to being distributed neatly in each DFT bin without double counting - and we can then use Parseval's Theorem to measure the total power (or the power across a band of frequencies) by summing the power in each bin. Windowing increases the equivalent noise bandwidth of each DFT bin (as given by the kernel, or discrete-time Fourier Transform, of the window), so that if we added up the power in each bin double counting would result and the total power would be over-estimated. The window itself also removes power from the signal since we are decreasing the signal over much of the time interval, resulting in a net loss (the loss from the window can be predicted based on the sum-square of the window samples, and the net loss is this loss plus the gain from the increase in kernel bandwidth). A single tone that is not spread out over multiple bins would only be effected by the power loss but not the gain from double counting, so the SNR is effected. This is called the processing gain of the window as detailed by fred harris in his classic paper "On The Use of Windowing in the Fourier Transform" and further detailed in this other posting:

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


Whats the technical term for a filter that allows a non unity stride when sliding it's window?

That's a decimating filter. As such, it serves as a 1/2-rate resampler.


Preprocessing audio signal for neural network classification

Some tips that might help:

Your FFT should be shorter than the features you are trying to detect.  An FFT will be calculated for a specific number of samples.  This represents a certain length of time (1024 samples at 1kHz sampling rate is 1024 milliseconds.)  If the time period of your FFT is longer than your signal features then they will all be smeared together - you won't be able to identify the features.  You will have to find a compromise for your FFT length between enough frequency resolution and enough time resolution.  The needed frequency resolution depends on the frequency content of features you need to detect, time resolution depndes on the length of the features (as well as the length of time when there are no features - which is itself perhaps a feature you need to take note of.)
Classify the features as to frequency range.  Make an average of the FFT-points for those ranges - this will give you the intensity of that feature.
Plot the intensity of the features against time.
Locate groupings of feature intensities that are unique for each type of heart defect - these are the things you need to be able to recognize.  For example: Pulmonic stenosis has low frequency (normal heart beat sound) with high frequency (murmur) at the same time, but no murmur type sound between beats - but it also has a click that you will need to detect.
You will need to detect the S1 and S2 by the pause between beats - long pause followed by short means that the last beat is the S2.
The examples on the site you referenced have been cleaned up, but they still have things that might try to sneak in as features.  I hear a 60Hz power hum, and some ringing from the filters that were used to clean things up.  There is also some sone background "white" noise in some of them.

I would try loading the files into Baudline and seeing how to identify the features and their combinations by eye first.  Baudline can show you how the freqeuncy content varies with time so that you can see the relationships between features and relative timing.  Use Audacity to convert the mp3 files to wav and to downsample to 1000Hz sampling rate before loading in Baudline.

The Baudline Color Aperture is handy for making certain parts of the diagram standout - this can make it easier to identify features.

Baudline is for Unix type systems only (Linux, BSD, Mac OSX) but I'm sure there's similar software for Windows.  Audacity is available for Windows as well as Linux and OSX.


Get spectral picture from a wavelet transform

Try this one:

m=cwt(y,1:100,'sym2','plot');
colormap(pink)

Result (btw, I cannot reproduce your curve in time domain with the function generate1):

Low scale values compress the wavelet and correlate better with high frequencies. While high scale values stretch the wavelet and correlate better with the low frequency content of the signal.


Can't get consistent magnitude response in MATLAB

As already mentioned by @Fat32 and @Florent Ecochard what your observing is a resolution issue with the produced plots. By default (e.g. using freqz without additional input arguments) 8192 points are used to plot the filter response; using less points, as done in you method 2, results in the observed behavior.

See the example code and figures below for an illustration:

%% Define Notch Filter
%
Fs     = 500;
Fnotch =  50;
BW     =   1;
Apass  =   1;

[b, a] = iirnotch(Fnotch/(Fs/2), BW/(Fs/2), Apass);
Hd     = dfilt.df2(b, a);

%% Plot Response

% method 1
% -> FREQZ default settings
%
[H1,F1] = freqz(Hd);
F1      = F1 * (Fs/(2*pi));  % convert to Hz / account for frequency normalization

% method 2
% -> Use frequency axis of 2048 points.
%
fs      = 500;
f_range = linspace(0,fs/2,2048);

[H2,F2] = freqz(Hd,f_range,fs);

% method 3
% -> Use frequency axis of 8192 points (same number as used in example 1).
%
fs      = 500;
f_range = linspace(0,fs/2,8192);

[H3,F3] = freqz(Hd,f_range,fs);

figure;
plot(F1,20*log10(abs(H1)),...
     F2,20*log10(abs(H2)),...
     F3,20*log10(abs(H3)) ...
    )
xlabel('f / Hz')
ylabel('|Hd| / dB')
legend('method 1 (freqz default)','method 2 (2048 points)','method 3 (8192 points)')

This then results in the following plots:


Why Low Frequency called Approximation and High Frequency Detail?

Intuitively speaking, anything that is 'high frequency' is something that is 'rapidly changing in time'. Anything that is 'low frequency' is something that is 'slowly changing in time'.

If you think about it, any time you have 'detail' in a signal or image, it means that you have, very quick, rapid variations in time or space. This then becomes the 'detail' of your signal/image.

In contrast, if you just care about how something changes generally, you care about long-term trends - you dont care about how things change with every sample, you just care about how things change over a long time, or over long space. In this case, you just want to 'approximate' your signal.


DSP or signal/image/data processing jokes

An airplane is leaving Warsaw (the capital of Poland), and gets caught in a terrible winter storm. The plane rolls, pitches and yaws. The crew is expecting the plane to crash or break up at any time. But one DSP student suddenly stands up and asks some passengers on the right aisle to move across the plane. Passengers are reluctant but the student insists. “Hurry, our fate depends on your move.” Most passengers quickly leap to the other side, and suddenly the plane stops shaking. The pilot, who had a DSP degree too,  enjoys calming everybody:

  "Thank you for the Poles who have moved to the left-half of the plane for stability."

Additional references:

The Routh Stability Criterion
Warsaw, Poland
Stability Revisited

Another visual pun (dedicated to @MarcusMuller) is this chocolet decomposition:


Can you present the convolution of sinusoidal with itself?


Removing drift from integration of accelerometer data

Double integration amplifies any offsets, non-linearities and noise.  These can't be removed without the use of some type of external reference point measurements (e.g. not from the accelerometer) or other information.  Calibration might be able to remove most of the ones that do not vary with time (temperature, etc.)

For instance, if the size of the room is known, you can robustly clamp position to be inside the walls and above the floor.


Implementation of Block Orthogonal Matching Pursuit (BOMP) Algorithm

The Block Orthogonal Matching Pursuit (BOMP) Algorithm is basically the Orthogonal Matching Pursuit (OMP) Algorithm with single major difference - Instead of selecting single index which maximizes the correlation we chose a set of indices, sub set of columns of the matrix and the solution vector.

A good reference for the algorithm is given in:

An Optimal Condition for the Block Orthogonal Matching Pursuit Algorithm.
Block Sparsity: Coherence and Efficient Recovery.

The code is given by:

function [ vX ] = SolveLsL0Bomp( mA, vB, numBlocks, paramK, tolVal )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveLsL0Omp( mA, vB, paramK, tolVal )
% Minimizes Least Squares of Linear System with L0 Constraint Using
% Block Orthogonal Matching Pursuit (OMP) Method.
% \arg \min_{x} {\left\| A x - b \right\|}_{2}^{2} subject to {\left\| x
% \right\|}_{2, 0} \leq K
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix (Fat Matrix). Assumed to be
%                           normlaized. Namely norm(mA(:, ii)) = 1 for any
%                           ii.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numBlocks         -   Number of Blocks.
%                           The number of blocks in the problem structure.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - paramK            -   Parameter K.
%                           The L0 constraint parameter. Basically the
%                           maximal number of active blocks in the
%                           solution.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - tolVal            -   Tolerance Value.
%                           Tolerance value for equality of the Linear
%                           System.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range [0, inf).
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  An Optimal Condition for the Block Orthogonal Matching Pursuit
%       Algorithm - https://ieeexplore.ieee.org/document/8404118.
%   2.  Block Sparsity: Coherence and Efficient Recovery - https://ieeexplore.ieee.org/document/4960226.
% Remarks:
%   1.  The algorithm assumes 'mA' is normalized (Each column).
%   2.  The number of columns in matrix 'mA' must be an integer
%       multiplication of the number of blocks.
%   3.  For 'numBlocks = numColumns' (Equivalent of 'numElmBlock = 1') the
%       algorithm becomes the classic OMP.
% Known Issues:
%   1.  A
% TODO:
%   1.  Pre Process 'mA' by normalizing its columns.
% Release Notes:
%   -   1.0.000     19/08/2019
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

numRows = size(mA, 1);
numCols = size(mA, 2);

numElmBlock = numCols / numBlocks;
if(round(numElmBlock) ~= numElmBlock)
    error('Number of Blocks Doesn''t Match Size of Arrays');
end

vActiveIdx      = false([numCols, 1]);
vR              = vB;
vX              = zeros([numCols, 1]);
activeBlckIdx   = [];

for ii = 1:paramK

    maxCorr         = 0;

    for jj = 1:numBlocks
        vBlockIdx = (((jj - 1) * numElmBlock) + 1):(jj * numElmBlock);

        currCorr = abs(mA(:, vBlockIdx).' * vR);
        if(currCorr > maxCorr)
            activeBlckIdx = jj;
            maxCorr = currCorr;
        end
    end

    vBlockIdx = (((activeBlckIdx - 1) * numElmBlock) + 1):(activeBlckIdx * numElmBlock);
    vActiveIdx(vBlockIdx) = true();

    vX(vActiveIdx) = mA(:, vActiveIdx) \ vB;
    vR = vB - (mA(:, vActiveIdx) * vX(vActiveIdx));

    resNorm = norm(vR);

    if(resNorm < tolVal)
        break;
    end

end

end

The MATLAB code is available at my StackExchange Signal Processing Q60197 GitHub Repository.
In the full code I compare the Block implementation to OMP to verify the implementation.


If you consider $X(s)$ as the only input of the system then the transfer function only exists for $T(s)=0$, and it is given by $G(s)$. Strictly speaking, with $T(s)\neq 0$, the system cannot be characterized by a transfer function because it is no longer linear. Note that a linear system must produce an output $Y(s)=0$ for $X(s)=0$. However, this is not the case for $T(s)\neq 0$.

You could also treat that system as a (linear) multiple input single output (MISO) system with two inputs $X(s)$ and $T(s)$. In that case there are two transfer functions:
$$H_X(s)=\frac{Y(s)}{X(s)}\Big|_{T(s)=0}=G(s)$$

and $$H_T(s)=\frac{Y(s)}{T(s)}\Big|_{X(s)=0}=3$$


In the code below I have manually created chirps where the frequency does not increase linearly, but behaves in a cosine-shape, as suggested by Atul Ingle . This waveform resembles the shape of the lines in the spectrogram nicely.

fs = 5000;
x0 = 1;
t = -x0:1/fs:x0;
yi = cos(t).*sinc(t/1.5).*cos(2*pi*500*(1.8-cos(2.2*t)).*t);
yj = cos(t).*sinc(t/1.5).*cos(2*pi*500*(0.85-0.5*cos(2.3*t)).*t);
subplot(2,2,1)
plot(t,yi+yj)
subplot(2,2,2)
plot(t,yj)
subplot(2,2,3)
spectrogram(yi+yj,512,511,256,fs,'yaxis')
subplot(2,2,4)
plot(t,yi)


How to know if the noise was removed in an image?

You can measure that with SNR, yes.  What you call 'signal' can get a little fluffy, but you should at least be able to compare two versions of the same image (before and after).

If your image has a nice homogeneous area in it, you can measure the variance of the pixel values in this region, and compare this before and after - so long as the image doesn't otherwise change.  Do make sure that you account for any scaling or shifting that a filter might be doing (an adaptive median filter shouldn't do either of these if it is implemented properly).

A common method for testing denoising filters is to start with a pure image (one that you define yourself), then add noise to it artificially then compare the denoised version to the original.  You will need to make sure to add the correct type of noise for your context (is it simple Gaussian, or Poisson, or ...?).


Difference between optical flow field and motion field?

well, have you ever seen a rotating barber's poll?

It looks like the stripes are moving up (optical flow), but of course the motion of the thing is a rotation in the horizontal plane.

The rationale really is that the optical flow is based on how the image changes. And if you just compared two frames from the animation above, you'd clearly say that the blue stripe moved upwards. However, if you'd open the device and looked inside, you'd see the motor that causes the mechanical movement (rotation).


You have the following problems:

First your $a$ and $b$ seems to be mis-named; replace b with a and a with b. Where $b$ is the numerator coefficients of the filter and $a$ is the denominator coefficients. Also note that the filter command works like this y = filter(b,a,x)...

Then also your $f_c$ parameter is not (directly) controlling the $-3$ dB cutoff frequency (in Hz) of the high pass filter that you intend to obtain. But it's controlling th eradian per sample cutoff frequency. So you bettwe use teh notation wc for that instaed of fc.

In any case, your input signal frequencies are so low that they are killed by the filter and hence you see no output. Just increase their frequencies to see them at the output. Incidentally you can also decrease $f_c$ well below $0.3$ to get a tighter highpass (a DC-notch) filter..

Based on your comment, since you want a cutoff frequency in radians/sample $\omega_c$  instaed of a analog Hz frequency, then the parameter $f_c$ yields the correct discrete-time cutoff frequency. You can view the frequency response of your filter by the line figure,freqz(b,a) . So the only problem remains is that your signakl frequencies ar too low. replace it like so

x1 = cos(0.1*fs*pi*n);
x2 = cos(0.3*fs*pi*n);

note that in your code the index $n$ has a scale of $T_s$ so the frequencies should be scaled by $f_s$ to get normalized frqeuencies in the discrete-time frequency.


Why do I have a high peak at the beginning of my FFT (not DC)?

You cannot simultaneously have a non-zero value at the first bin of the DFT (the 0 component) and a zero-mean since the first component is by definition of the DFT the (scaled) mean of the signal:
from the DFT:
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}$$
$$X[0] = \sum_{n=0}^{N-1}x[n]e^{j0}= \sum_{n=0}^{N-1}x[n] = N\mu_x$$
So there is a mean in the digital data, and perhaps the OP's confusion is there is not a mean of any significance in the analog waveform being sampled. There are several opportunities to introduce offsets notably in the ADC converter itself, or it may be as simple as the mapped amplitude range of the digital waveform includes an offset. But taking the mean of the data that is being FFT'd should be consistent with what is above, with further scaling by any windowing that may be done according to $\sum_{n=0}^{N-1}w_n$.


How to determine the number of multipath?

You can't really tell that from anything. The maximum delay spread together with your sampling rate, however, gives you the number of paths that you can tell apart with a single observation of the channel. If you model every sample of your discrete-time channel impulse response as a path.
That's not really a sensible "channel model", however: the tapped delay line model works well enough for point reflectors and static channels. However, your channel changes, and your reflectors are moving: you get paths as "smeared" peaks in your ambiguity function.
Still, the tapped delay line model with each sample of the channel impulse response being one tap is immensely useful, as it allows you to equalize the channel by OFDM.
Many scenarious consider point reflectors (which sadly don't always sit "exactly" sampling periods apart) that are moving, thus "smearing" and leaking across multiple taps. As soon as you add anything that starts to track different path components, you might be able to sense more paths than you have taps!
So, honestly, there's no "maximum number of paths" that you can distinguish if you have enough time and if they are uncorrelated enough.


Should I sample at twice the bandwidth or twice the highest frequency?

Well, according to Nyquist-Shannon theorem you should sample at a rate which is at least twice the highest frequency you want to capture. This is also referred to as the sampling theorem because it "forms the basis" for sampling.

Now regarding your specific question, I have to say that cos(4*pi) (along with the rest of the components) is just a number which is equivalent to a DC signal. In case you imply cos(2 * pi * 4 * n/fs) (the same for the rest of the components), where n is the sample index (constantly increasing) and fs is the sampling frequency, you should sample at twice the highest frequency (which in this case is 4), which is 8.

As you might be able to say, this would produce a result (just for the first component, but similar for the rest) cos(n * pi), which has a period of 2..., which in turn means that you would get precisely two samples per period. Now in the case you happen to sample at the zero cross, all samples would end up being zero (just for this highest frequency).

I strongly suggest you have a look at some academic literature regarding sampling. I would suggest Proakis' & Manolakis' Digital Signal Processing - Principles, Algorithms and Applications and Orfanidis' Introduction to Signal Processing


How to calculate pitch fundamental frequency f( 0) ) in time domain?

This paper by Philip McLeod details finding pitch in the time domain, and specifically that it avoids the reduced resolution of equivalent frequency domain approaches due to the windowing required to eliminate spectral leakage.

https://cdn.instructables.com/ORIG/F3J/3Q7U/J20P2R5T/F3J3Q7UJ20P2R5T.pdf

I don't have experience with this approach but would be very interested in comparative results to the best frequency domain approach- if you could post them, assuming you pursue this further.


Wavelet Transfrom + Power Spectral Density (using Matlab)

You need to do 1-Dimensional wavelet package decomposition. Here is the matlab documentation for the command to do this.

Here is some example code based on what you've described above:

wpt = wpdec(X, 2, 'db1', 'shannon')
plot(wpt)

This will provide a nice clickable graph. As you'll see in the documentation, you can provide a lot of different parameters for the decomposition.


White noise gain of uniform spline interpolation

This figure is the squared norm or energy of the impulse response:

$$
\|h(t)\|^2=\int_{\mathrm{R}} |h(t)|^2 dt
$$

As you point out correctly, this is the noise amplification. Considering a signal

$$
r(t) = x(t) + n(t)
$$
with $x(t)$ is the signal of interest and $n(t)$ is AWGN with variance $\sigma^2$. Then, after filtering

$$
\begin{align}
y(t) &= h(t)*r(t)\\
&=h(t)*x(t) + \int_{\mathrm{R}} h(t-\tau)n(\tau)d\tau\\
&=h(t)*x(t) + \tilde{n}(t)
\end{align}
$$

with $\tilde{n}(t)$ being the filtered noise. Now, considering

$$
\begin{align}
E[\tilde{n}(t)\tilde{n}^*(t)] &= \int_\mathrm{R}\int_\mathrm{R} h(t-\tau_1)h^*(t-\tau_2)E[n(\tau_1)n(\tau_2)]d\tau_1 d\tau_2\\
&=\sigma^2\int_{\mathrm{R}^2}h(t-\tau_1)h^*(t-\tau_2)\delta(\tau_1-\tau_2)d\tau_1d\tau_2\\
&=\sigma^2\int_\mathrm{R}h(t-\tau_1)h^*(t-\tau_1)d\tau_1\\
&=\sigma^2\|h(t)\|^2
\end{align}
$$

where we have used that the noise is white, i.e. its autocorrelation is a delta-function. Note that after filtering generally the noise is not white anymore, but colored.


Sobel and Prewitt operators decomposition

The part:
\begin{bmatrix}
1\\
1\\
1
\end{bmatrix}
acts on vertical pixels
\begin{bmatrix}
x_{1,\cdot}\\
x_{2,\cdot}\\
x_{3,\cdot}
\end{bmatrix}
as a weighted sum with equal weights, and if divided by $3$ like a 3-point average $(x_{1,\cdot}+x_{2,\cdot}+x_{3,\cdot})/3$.

Then,
\begin{bmatrix}
+1 & 0 & -1
\end{bmatrix}
acts on a row of horizontal pixels
\begin{bmatrix}
x_{\cdot,1}&
x_{\cdot,2}&
x_{\cdot,3}
\end{bmatrix}
like the discrete 3-point derivative: $(x_{\cdot,1}-x_{\cdot,3})/2$ up to the $2$ factor. Here, up to a global $6$ factor, it acts like a smoother vertically combined with a gradient horizontally. The same reasoning works when you switch directions. This notion relates to the separability of a 2D mask into the product of two 1D operators.

Additional sources:

What's logic behind the constuction of Sobel's filter in image processing?
How to find out if a transform matrix is separable?