Estimate Instantaneous Frequency Using LMS Algorithm I hope someone can help me with the following problem:

I want to estimate the frequency of a sound file that is composed of a sinusoidal with varying frequency and additive white noise:

$$ x \left[ n \right] = \sin \left( \varTheta \left[ n \right] \cdot n \right) + w \left[ n \right] $$

Now I implemented a linear predictor using LMS algorithm to estimate the parameters of an AR-Process (As described in Haykin: Adaptive Filter Theory).

$$ x(n) = -\sum_{k=1}^M a_k(n)x(n-k) + v(n) $$
$$\hat{w}(n+1) = \hat{w}(n) + \mu x(n-k)f_M(n) $$

With $f_M(n) = u(n) -\sum_{k=1}^M \hat{w}(n)x(n-k)$ being the prediction error.

From that I calculate the time-varying frequency function:
$F(\omega,n) = \frac{1}{|1-\sum_{k=1}^M \hat{w}_k(n) e^{-j\omega k}|^2}$

This looks like the following:

For now I compute the instantaneous frequency by searching for the maximum of F for each time sample n. But the problem is, that it is quite noisy and the frequency overshoots a lot at the transitions:

Now I wanted to ask if there are any ways to optimize the computations, I would be very happy for some ideas.

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

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

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

where

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


Spurious When $F_s/F$ is not an integer We are using a 16-bit DAC for a waveform generation between  to $500\textrm{ MHz}$ with the sampling frequency of $1200\textrm{ MHz}$. The specification for the waveform generation is $10\textrm{ Hz}$ frequency resolution throughout the band  $>60\textrm{ dB}$ SFDR. The spectral performance of the DAC is fine (more than $65\textrm{ dB}$ SFDR) for the frequencies which is satisfying the condition:

$$\frac{F_s}{F}=k \quad\big\vert\quad k\in \mathbb Z$$

The frequencies which is not satisfying the above condition yields some spurious which is not meeting my SFDR specification. Since the samples of every cycle will not be same, it's creating an envelope over the waveform, due to a phase values repeating itself periodically over a number of cycle. This envelope causing a spurious.

If we generate a frequency of $200.01\textrm{ MHz}$ in MATLAB with $F_s = 1.2\textrm{ GHz }$, we can clearly notice the envelope with periodic time period of $16.65\textrm{ $\mu$s}$ ($60\textrm{ kHz}$). The figure below is the MATLAB time-domain plot for $200.01\textrm{ MHz}$.

With my DAC, I'm getting a spurious at $\pm 60\textrm{ kHz}$ and $\pm 30\textrm{ kHz}$ apart from my generated frequency($F$).

The image below is the DAC output spectrum for $200.01\textrm{ MHz}$ captured with spectrum analyzer.

Is the spurious observed is due to the envelope?

MATLAB Simulation Code;

Fs = 1200e6;    %Sampling Frequency
T = 1/Fs;       % Clock Period
N = 2^18;       % No. of samples
t = (0:N-1)*T;  % Linear time interval
F = 200.01e6;   % Required Frequency
phase = 2*pi*F*t; % Consecutive Phase values in rad
sig = sin(phase); % Time Domain Signal
figure;
plot(t, sig), xlabel('Time in Sec'), ylabel('Amp'), title('Time Domain Plot for 200.01MHz');
wrap_phase = mod(phase, 2*pi);
figure;
plot(t,wrap_phase), xlabel('Time in Sec'), ylabel('Phase in rad'),           title('Phase values for 200.01MHz');

Suggestions are greatly appreciated.

Loganathan N

As people who responded in DSPrelated to your question, I believe the fluctuation which are seen are MATLAB plot artifacts not a real envelope. Since you doubt if Fourier analysis could help us to see the spurious signal, I perform a more elaborate spectrum sensing method (MUSIC) to identify if there is other frequency contents in the signal through the code below:

s = 1200*10^6;%Sampling Frequency
T = 1/Fs; % Clock Period
N = 2^12; % No. of samples
t = (0:N-1)*T; % Linear time interval
F = 200.01*10^6; % Required Frequency
phase = 2*pi*F*t; % Consecutive Phase values in rad
sig = sin(phase); % Time Domain Signal
%% SPECTRUM ESTIMATION
X=corrmtx(sig,200,'mod');    % MUSIC SPECTRUM ESTIMATION
[W,P] = rootmusic(X,4);
f=Fs*W./(2*pi); %Frequencies
unique(abs(f)) %Show Frequencies
P              %Power of those frequencies

if you check out the P for each frequency component, you'll see only the main frequency (here 200.01MHz) has non-zero value.


Feature vector for MFCCs I am working on a project, for which I am finding the onset times of an audio and then extracting the MFCC coefficients for all the times.

Thus if I have n onset times, I would have n feature vectors, corresponding to the MFCC coefficients at those points.

Now I want to characterize a "bol" of the tabla by the the waveform in between the start of one onset and the start of the next onset. This would essentially mean that for one onset time, I would have to take the MFCCs feature vectors at multiple points, and I would have a matrix of onset values for each onset point, instead of just a single vector at each point. Also, the matrix wouldnot always be of the same size, since the gap between 2 onsets might not always be the same.

How do I combine these multiple MFCC vectors to one, or how to I use this matrix for machine learning?

If I understand the problem well, you have a segment of audio between two onsets, for example between an onset at t=2.0s and an onset at t=2.8s, and your MFCC function returns many MFCC vectors - corresponding to small analysis windows between these two instants (In our example, 80 vectors if your MFCC function is set for a rate of 100 analysis frames per second).

But you want a single vector "summarizing" the MFCC for the entire region between the two onsets.

Option 1:

Increase the window size / decrease the frame rate of your MFCC routine so that it yields fewer vectors.

Option 2:

Build a feature vector with the mean and standard deviation of all the MFCC vectors collected within your time segment.

In both cases, set an upper bound on the size of the segment analyzed by such a process - I'd say at most 200ms. For example, if you have an onset at t=5.0s and the next onset is at t=7.0s, the percussive event is likely to be located within the first couple of hundred ms (say between t=5.0s and t=5.5s) - the rest being silence or background noise.


What’s a good book on multi-rate dsp? Was Studying a dsp design recently with many decimation and interpolation filters, and was wondering if there is a recommended book specifically on the subject of multi-rate dsp?  Preferably a book with worked examples and exercises at the end of the chapter.

Assuming you have completed the basic course in DSP with a book like Oppenheim's Discrete-Time Signal Processing book then the next step could be anyone of the followings:

Multiresolution Signal Decomposition: Transforms, Subbands,Wavelets, A.Akansu

Multirate Digital Signal Processing, Rabiner

Multirate Systems and Filter Banks, P.P.Vaidyanathan

Wavelets and Filter Banks, Gilbert Strang

Multidimensional Digital Signal Processing, Dudgeon

Subband Adaptive Filtering: Theory and Implementation", by Woon-Seng Gan, Kong-Aik Lee, Sen M. Kuo

Wavelets and subband coding, Vetterli

Advanced DSP, John G Proakis

For a good book on multi rate DSP from communications perspective, please refer to the book below

Multi rate signal processing for communication systems, Friedrich J Harris

For a MATLAB specific handling, please refer to the book below

Multirate Filtering for Digital Signal Processing: MATLAB Applications, Ljiljana Milic

Personal recommendation the books by Vetterli and Fredrich Harris


Estimating device height from 3-axis accelerometer data (possibly by supervised learning) I'm trying to estimate the height of my device using an cheap integrated Freescale accelerometer.

First, I estimate the earth direction when measured acceleration is close to 1g, then I project the measurement on that direction and integrate twice.
Very quickly the estimated height drifts.

I was thinking maybe to use some machine learning algorithm, by feeding it aligned height from video tracking as examples.

Can you suggest algorithms that would help doing that?
Or any other way to make height estimation more precise?

Edit: Motion is known to be limited in a range of ±0.5m.

Any tiny offset error will cause a double integration to fly off quadratically in time.  Even if you manage to calibrate out this error for one sample test setup, any drift in offset due to change in temperature, orientation, age, power fluctuations, etc., will eventually render that calibration useless.

Continuous recalibration (airspeed, compass, optical fix, etc.) is likely required for a more useful result over a sufficiently long time period.


dB SPL readings differ from external sound level meter? I'm trying to write a program that takes values from a MEMS microphone and calculates the dB of the incoming audio data.

The microphone i'm using is this: https://www.mouser.es/datasheet/2/218/ph0645lm4h-b-datasheet-rev-c-1525723.pdf using the I2S protocol, 24 bits signed integers.

The problem i have is that my dB SPL calculations are reading at 80 dB SPL at silence whereas my sound level meter reads around 30/40 dB(A).

Here is how i'm calculating the db SPL of the microphone input:

#define DC_OFFSET 419430

// this int32_t contains a 24 bit signed integer
int32_t sample = GetSample();

// Remove the DC offset.
sample = sample + DC_OFFSET;

sample = sqrt(pow(sample, 2));
double sampleNormalized = ((double)(sample) / 8388607);

// Convert the normalized samples to dB FS
double dbFS = 20.0 * log10(sampleNormalized);
double dbSPL = dbFS + 120.0;

Why do my readings read from 80 dB SPL whereas my sound level meter reads from 30/40 dB(A)?

UPDATE

Now my GetSample() function returns a signed 18 bit values, a few test samples results in these values, before DB offset adjustment:

78528
77376
73472
80384
79872
83648

DC offset defined as (+-100 5%), max value of 18 bits is 131071, therefore DC offset is (131071 * 0.05) = 6553.0.

When i pass it through this algorithm to convert it to dBSPL, the values now come out even higher than before, but the value responds how it should.

double ref = 131071.0;
double dc_offset = ref * 0.05;

// read 18 bit signed value
int32_t sample = GetSample(); // 18 bits wide.

// remove the dc offset.
double adjustedSample = sample - dc_offset;
adjustedSample = fabs(adjustedSample);

// normalize. the sample to the range 0.0 to 1.0.
double sampleNormalized = ((double)(adjustedSample) / ref);

// map the range from 0.0->1.0 to -1.0->1.0.
sampleNormalized -= 0.5;

sampleNormalized = abs(sampleNormalized);

// convert the normalized samples to dB FS
double dbFS = 20.0 * log10(sampleNormalized);
double dbSPL = dbFS + 120.0;

I'm still calculating db SPL values that are 88 dB SPL at quiet background, but the value actually responds correctly, i will apply the A filter as @Max posted and return the results afterwards.

Your calculations are correct, although you should do some kind of sliding average over a few hundred samples or so. The problem seems to be that the sound level meter employs an A-weighting, while your code does not. If there is low frequency noise in your measurement environment, the readings at "silence" can easily differ by 40-50dB. Try measuring a test tone at 1kHz, results should be closer together. Also, have a look at

https://en.wikipedia.org/wiki/A-weighting

You can see the huge low frequency differences there.


pole/zero locations for real and imaginary signal In Z-Transform,

For a real signal,
 $x(n)$ =$x^*(n)$ .
Taking Z-transform on both sides,

$X(z)$=$X^*(z^*)$ ,  which gives certain pole/zero condition

similarly for a purely imaginary signal            ,

$x(n)$ =$-x^*(n)$ , if we take Z-Transform on both sides , $X(z)$=$-X^*(z^*)$

which gives similar pole/zero condition as that of real signal(poles and zeros are not effected by negative sign).

But  My professor says if the poles/zeros are symmetric with respect to real axis, then it is a real signal. But i am getting same condition for both real and imaginary signal.

Can anyone help me? (examples might be helpful)

I think there is a confusion between signals and systems here. A signal doesn't have any poles and zeros. Only systems have poles and zeros.

A system is a thing that turns one signal into another one. For example a loudspeaker turns a voltage signal into a sound signal. The voltage going in is a signal, the sound coming out is a signal, but the loudspeaker itself is a system.

Linear Time Invariant Systems can be characterized by their transfer function or their impulse response. The transfer function in the Z domain can be written as a rationale function, which in turn can be characterized through poles and zeros.

If the system has only a real valued input and only a real valued output, the impulse response must be real as well and that in turn will require the poles and zeros to be either real or be conjugate complex pairs.


Simple questions related to dft I'm reading a book about dsp and doing some exercises.
Here is question that confused me:

A peak appears at index number 19  when a 256 point DFT is taken of a signal

1) What is the frequency of the peak expressed as a fraction of the sampling rate? Do we need to know the actual sampling rate to answer this question?

2) What is the sampling rate if the peak corresponds to 21.5 kHz in analog signal?

My answers:

1) 19/256; no we don't need to know actual sampling rate

2) I think that we can't answer this question without additional info. Or maybe it is  ((21.5/19)*256)*2 kHz?

Thanks in advance

The actual frequency in Hertz corresponding to DFT index $i$ with a DFT length of $N$ is

$$f=\frac{i}{N}f_s,\quad i=0,1,\ldots N-1\tag{1}$$

where $f_s$ is the sampling frequency. So if indices are numbered starting from $0$ (and index number $19$ is indeed the $20^{th}$ frequency bin, as assumed in (1)), then your answer is correct.

As for part 2 of the question, the answer is also given by (1), so

$$f_s=\frac{Nf}{i}$$

which gives $f_s=289.7\text{ kHz}$ for the given values of $i$, $N$, and $f$.


How transmission speed and bandwidth are linked? I don't understand why if I have a larger bandwidth I can transmit data faster.
Is this linked with this property of fourier transform?

thanks.

Actually you also have to consider the noise floor of the communication channel and transmitted signal power to define the limit of data rate through an AWGN channel. The information capacity $C$ (bps) of the channel is given by the Shannon's formula:

$$ C = B \cdot \log_2( 1 + \frac{ \sigma_x^2}{ \sigma_n^2 } ) $$

where $B$ is the channel bandwidth in Hz., and $\frac{ \sigma_x^2}{ \sigma_n^2 }$ is the channel SNR (Signal to Noise Power Ratio).

One ideal case happens when there's no noise, then the capacity goes to infinity. This means that on a noiseless (ideal) channel you would transmit all the information (and more) instantly.


calculating estimate of a state for a system with two observations of the state from different times I'm struggling with a problem that I just can't seem to get a grasp of. I'm supposed to calculate an estimate for state $x(k)$ at times $k=1,2,3$ from the state-space

$
\begin{align}
x(k+1) &= A  x(k) + B  u(k) + w(k) \\
y_1(k) &= x(k) + v1_(k) \\
y_2(k) &= x(k-1) + v_2(k) \\
\end{align}
$

where the covariances of the noises are given as

$E[w(k)w(k)']$ , $E[v_1(k)v_1(k)']$, $E[v_2(k)v_2(k)']$

and

$u(k)$ is given for $k=1,2,3$

$y_1(k)$ is given for $k=1,2,3$

$y_2(k)$ is given for $k=2,3$

What kind of technique am I supposed to use? Information Kalman filter?

HINT

$$ \left| \begin{array}{cc} A & 0 \\ 0 & 1 \end{array} \right| $$


Is $\cos(\pi \sqrt n)u[n]$ stable? I took a z transform and got a double pole at $z=1$, but I don't know if that's correct.

I'm  lost because I don't know if $\cos(\theta)$ converges or diverges or what that means for $h[n]$ being absolutely summable.

The system with impulse response given by $h[n] = \cos(\pi\sqrt{n})u[n]$ is BIBO-unstable because the sum $\sum_{n=-\infty}^\infty |h[n]]$ diverges instead of being convergent as is needed for BIBO-stability. Note that
for all positive integers $k$, $h[k^2]=\cos(\pi k)$ has value $\pm 1$ and so
$$\sum_{n=-\infty}^\infty |h[n]] = \sum_{n=0}^\infty |\cos((\pi\sqrt{n})|$$
is a sum that contains infinitely many $+1$ terms (and all the other terms are guaranteed to be positive too since $\cos(\pi r)=0$ if and only if $r = k+\frac 12$ where $k$ is an integer, and there is no integer $n$ whose square root is of the form $k+\frac 12$). So, $\sum_{n=-\infty}^\infty |h[n]]$ diverges, and the system is BIBO-unstable.


Understanding the Gabor filter function I need to implement a script for generating features from an input image by using the Gabor filter. I have no past experience of wavelets and I'm just learning Fourier analysis (I understand the basic idea behind Fourier analysis and transform) so they can't help me to understand Gabor filter, because I need to have an implementation done in a week. I don't need to understand the foundations of the Gabor filter function, but I would want to understand to some extent of what it is and what does it do. What are the parameters? What do they mean? What is the output of the function? For example this is the formula I copied from Wikipedia:

$$g(x,y;\lambda, \theta, \psi, \sigma,\gamma) = \exp\left(-\frac{x'^2+\gamma^2 y'^2}{2\sigma^2}\right)\exp\left(i\left(2\pi\frac{x'}{\lambda}+\psi\right)\right)$$

Now my obvious question is: What does this mean? What does the variables mean? According to Wikipedia:

$x, y$: I assume these coordinates specify the pixel value of an image at coordinates $(x,y)$ (2. This is OK, understood)

$\lambda$: represents the wavelength of the sinusoidal factor (Sinusoidal factor, huh? 3a. How do you select it? 3b. Where does it come from? 3c. Is it an arbitrary number or what? Freely chosen?)

$\theta$: represents the orientation of the normal to the parallel stripes of a Gabor function (4. What does that mean?)

$\psi$: is the phase offset (5. offset of what? How is this value determined? Is it freely chosen?)

$\sigma$: is the sigma/standard deviation of the Gaussian envelope (6. Need more explanation...)

$\gamma$: is the spatial aspect ratio, and specifies the ellipticity of the support of the Gabor function (7. Again need more details and more explanation)

And most importantly:

$$g(x,y;\lambda, \theta, \psi, \sigma,\gamma) = X$$

What is the output value $X$? What does it mean?

As I mentioned I don't need thorough explanation of the theory, because I bet it is long and reading a 1000 page book on unknown subject is not an option for me right now. I need to have have a black-box understanding of this function so that I can implement it in code and most importantly understand what is the input and what is the output.

Thank you for any help!! =)

P.S.

I read this post:

https://math.stackexchange.com/questions/259877/value-of-x-y-in-computing-gabor-filter-function

but it doesn't answer my question thorough enough :)

A matlab implementation is in this answer: https://dsp.stackexchange.com/a/14201/5737

1) The wikipedia formula is a little bit too general.

2) If you know the basics of the with Fourier transform then you empirically know that: The image is viewed as being formed by superimposing a series of sinusoidal waves of various frequencies oriented in all kind of directions. Each "pixel" in the transform tells us the "intensity" of such a wave. The position of the "pixel" tells us the frequency and orientation of the wave. In practice, one wants to select only certain waves, having a specific frequency and a specific orientation.

So there you have it: The Gabor transform is one of many so called band pass filters that allows you to "cut" the Fourier transform and isolate only specific information. Another important information is that each Fourier "pixel" is a complex value (real and imaginary part)

3) Parameters: Two parameters are already shown:

3.a) The tuning frequency $f_0$, or tuning period $P_0$, or $\lambda$ establish to what kind of sinus wave the filter will respond best. ( $f = 1/P_0 = 1/\lambda$ or $f=\pi/\lambda$ depending on the specific implementation ) Basically, a smaller $P_0$ means a denser sinus wave. A larger $P_0$ means larger waves. $P_0$  is in pixels (3, 5, 30, etc. pixels) Don't go under 3 pixels or beyond $W/2$ or you get nasty effects. $W$ is the width of the image, if the image is square. You specify this using P0 parameter in the matlab code.

3.b) Central angle. These waves can have any directions. You want to select only waves at specific angle. So, the second parameter is the tuning angle, $\theta_0$ or $\theta$ in your formula. Usually, in radians. Orient in matlab code.

One cannot isolate only a certain frequency or a certain orientation. (Search for uncertainty principle in the textbooks. Yes, is similar to the one in physics) But one can tune how much of the nearby frequencies will leak. The next two parameters specify that:

3.c) $\Delta F$, frequency bandwidth expressed in octaves. Usefull values, 1.5, 2, 3. Larger values means capturing a broader range of frequencies. There is a price for a tighter bandpass, a poorer spatial localization. Why, again, textbook. FBW in matlab code.

3.d) $\Delta \theta$ the angle bandwidth. Expressed in radians. $\pi/3$ or $\pi/2$ works just fine. ABW in matlab code.

The relation between $\Delta f$, $\Delta \theta$ from matlab code and $\sigma, \gamma$ from Wikipedia have a formula but is not essential for understanding Gabor. $\psi$ is again not important for basic understanding.

So you have it: $X$ from Wikipedia is a 2D matrix of numbers containing a convolution mask. You take the original image, filter it with the convolution mask and get another image. This new image is the "Gabor response" for the original image.

The matlab code does these two steps together. Constructs the Gabor filter with the specified parameters and performs convolution. The results  ReConv,ImConv are the responses. Each  response "pixel" have a real and an imaginary part. If you want to use this code you usually must compute the energy for the response: $E = \sqrt{a^2+b^2}$ where $a$ is the real part of the response (ReConv) and $b$ is the imaginary part (ImConv), for each pixel.

There you have it:
1) Build a Gabor filter specifying $P_0, \theta_0, \Delta f, \Delta \theta$
2) Convolute your image with the filter. You will get two values per each pixel.
3) Compute the energy $E$ and get the intensity of the response for each pixel in the original image

Another intuition: Suppose you want to select edges stretching on an orientation perpendicular to $\pi/6$ and a certain width of 20 pixels. You can build a Gabor filter with $P_0 \approx 20, \theta_0 = \pi/6, \Delta f = 2, \Delta \theta = pi/2$. TAKE CARE that there is no equal relation between the dimension of your edge and the $P_0$ parameter. You should try various values and see what works best. The other parameters ($\Delta f, \Delta \theta$) touch them when you have some experience tuning the first two.

Hope it helps!

Cristi

Update

Here is a site which allows you to play a bit with Gabor parameters and note the results: http://www.cogsci.nl/pages/gabor-generator.php

After a quick look, both freqnency and angle bandwidths are tied together as "Standard Deviation in pixels ... to the Gauss envelope". The rest of the parameters are easily identified. However, note that there are numerical differences! Eg. the phase can be expressed in 0-1 interval or 0-$\pi$ interval.


How to express FM and PM the most simple way? I'm trying to understand how to tell phase modulation from frequency modulation by looking at their respective expressions, simplified, particularly frequency modulation. Also, it would be nice for someone to illustrate mathematically the difference between vibrato and FM synthesis.

In the simplest terms, the difference between analog Phase Modulation (PM) and Frequency Modulation (FM) can be defined according to which of the two (phase or frequency) is being directly proportional to the message signal as outlined by the following:

An angle modulated signal is defined as: $$x(t) = A_c\cos(\omega_c t + \phi(t)) $$ Where $\theta(t) = \omega_c t + \phi(t) $ is the instantaneous angle  of the modulated cosine, $\phi(t)$ is the instantaneous phase and  $ \omega_i = \frac {d\theta(t)}{dt} $ is the instantaneous frequency

Now consider the phase being given by $\phi(t) = K_p m(t)$ , where $m(t)$ is the message signal and $K_p$ is a simple constant for phase modulation. As clearly evident, here the Instantaneous Phase, $\phi(t)$, is directly proportional to the amplitude of the message signal $m(t)$ and hence this type of angle modulation is called as Phase Modulation (PM)

On the other hand, consider a message signal $m(t)$ such that $\phi(t) = K_f \int {m(\tau)d\tau}$, where $K_f$ is the constant for frequency modulation, then the instantaneous frequency is: $$ \omega_i = \frac {d\theta(t)}{dt} = \omega_c + K_f m(t) $$
As can bee seen, now instead, the instantaneous frequency in radians is directly proportional to the amplitude of the message signal $m(t)$ and this type of modulation is, therefore, called as Frequency Modulation (FM)


Frequency tracking with huge noise I'm working on a frequency tracking problem with noise, where the amplitude of the noise is orders of magnitude higher than the amplitude of my signal (~1000x). Some details:

The signal is roughly sinusoidal, less than 1 Hz with non-stationary frequency and amplitude
The noise is caused by motion, and is (in many cases) not periodic
I have a 3-channel accelerometer reference for the noise. My target signal has some correlation with motion, so adaptive filtering may have issues

Empirical Mode Decomposition and even simple windowed FFT work well in the no-motion case, but fall apart as soon as motion kicks in - the broadband effects of the noise mask the frequency-domain peak of the signal. I tried RLS adaptive filtering, but I believe the aperiodicity of the noise prevents the filter from converging to anything useful. Would something like Extended Kalman Filtering be able to handle this case? What other options are there when your noise overpowers your signal?

Update

Phase-locked Loops:

According to this article, PLLs...

  ... can demodulate FM signals with very high accuracy and reliability, or it can detect signals buried in noise, but it can't do both in a single configuration, because the two tasks require very different setups and assumptions. In its role as an FM detector, a PLL doesn't reject noise very efficiently, and as a weak-signal detector, it can't decode FM very efficiently. The reason should be obvious — to detect FM modulation of a given bandwidth, the PLL's feedback loop low-pass filter (Figure 2, green) must be opened up enough to allow the modulation's bandwidth to pass unimpeded, but this causes the PLL to become more susceptible to noise.

Kalman Filtering:

I've been testing out an Unscented Kalman Filter (UKF) for the problem. I found examples of people using Extended Kalman filters, which work ok for the simplest problem of tracking a sine wave because it's easy to derive the Jacobian of a sinusoid. I opted to try the UKF because of the ease of incorporating further non-linear constraints / future problems where the Jacobian isn't easily accessible.

I set up the filter using a state variable of $\begin{bmatrix} A & \omega & \phi & b\end{bmatrix}^T$, where the measured signal associated with the state is $z = A \sin (\omega - \phi) + b$.

The process model assumes that $A$, $\omega$ and $b$ remain constant, and the phase advances according to the current frequency: $\phi[t+1] = \phi[t] + \omega[t] \times dt$ ($\phi$ is constrained to $[-\pi, \pi)$).

The problem I find with this model is that solutions are not unique. For instance, in the case pictured below, the filter has converged by sending the frequency to 0, and manipulating the phase variable to influence the signal:

In each subfigure, the black line represents the estimate from the filter. The top subfigure is the fit of the model against noisy measurements (blue dots). The remaining subfigures show each variable of the state vector - colored dots are the true value that we would like the filter to converge to.

Even if I could restrict the argument of the sin to one process variable, it would always be possible to set the frequency to 0 and use the amplitude or baseline shift to fit the signal.

My sense is there are clever ways of designing the process model to constrain the behavior of the filter. Can anyone suggest a better way to set up the problem?

For reference, I have attached a python 3 notebook with the code I am using.

might look at:

B. F. La Scala and R. R. Bitmead, "Design of an extended Kalman filter > frequency tracker," Acoustics, Speech, and Signal
Processing, 1994.

  Abstract: This paper discusses the problem of designing a frequency tracker using the extended Kalman filter (EKF). The design
  trade-off between balancing noise rejection and tracking at a maximal
  slew rate is considered, as are the effect of varying design noise
  covariance values. The performance penalties for over- and
  under-design of noise covariances are examined and theoretically
  supported design guidelines are suggested.

There's also the book :

Quinn, Barry G., and Edward James Hannan. The estimation and tracking of frequency. Vol. 9. Cambridge University Press, 2001.

As Peter noted, in SONAR, tracking low SNR tones is done but they typically run off a DFT.  If you want to go this route, you would have to do some experimentation on bin widths and short term averaging.  Fundamentally, you need to filter the signal down to a band that results in a positive SNR.  The band can contain more than a few bins, but the signal has to display a reliable peak when it occupies the bin.  A much more expensive approach would be a bank of chirps.

I would also suggest, for your remark about process and measure noise being correlated.  One of the reasons I like the book is that they derive the KF directly from the case where measurement and process noise is correlated.

Kailath, Thomas, Ali H. Sayed, and Babak Hassibi. Linear estimation. Vol. 1. Upper Saddle River, NJ: Prentice Hall, 2000.


Understanding Matlab's Automatic Gain Control function I am trying to understand the automatic gain control block provided by the communications toolbox. The documentation is here: Documentation

My questions are in regards to two of the configurable parameters, the
DesiredOutputPower and MaxPowerGain.

Given a value X for DesiredOutputPower, how does one compute what the reference value A is?

Given a value X for MaxPowerGain, how does one compute what the upper limit for g is?

For example:

If DesiredOutputPower = 2 then A = 0.693147180559945

If DesiredOutputPower = 6 then A = 1.791759469228055

And

If MaxPowerGain = 10 then the upper limit for g is 1.15129255

If MaxPowerGain = 60 then the upper limit for g is 6.90775527

The documentation in regard to DesiredOutputPower states:

  Specify the desired output power level as a real positive scalar. The
  power is measured in Watts referenced to 1 ohm. The default is 1.

And in regard to MaxPowerGain:

  Maximum power gain in decibels

  Specify the maximum gain of the AGC in decibels as a positive scalar.
  The default is 60.

How are these numbers determined?
What are the formulas? I've been trying to apply the formulas I read here:
The dB in Communications
but I havent been able to figure it out, please help.

Thank you!

From the diagram in the Algorithms section of the documentation you can see how the different quantities are computed:

Note that $z$ in the diagram is an estimate of the output power.$^1$ The error signal $e$ is computed by comparing the reference value $A$ to $\ln(z)$. So if you choose $$A=\ln(P)\tag{1}$$ then the average output power will be adjusted to the specified value $P$.

Concerning the maximum gain, note that the input is multiplied by $\tilde{g}=e^g$. If $G$ is the maximum gain in dB, you have the following relationship:

$$G=20\log_{10}(\tilde{g})\tag{2}$$

i.e.,

$$\tilde{g}=10^{G/20}\tag{3}$$

and, consequently,

$$g=\ln(\tilde{g})=G\cdot \frac{\ln(10)}{20}\approx G\cdot  0.11513\tag{4}$$

1. There is an error in the diagram: the output of the detector must be multiplied by the square of the gain (because it's an estimate of the power). This is correctly represented by the equations below the diagram in the documentation.


Drawing a Spectrogram off of FFT - what are the decibels relative to? I am trying to draw a spectrogram from FFT data I get from audio files.
The problem is that I can't figure out is what the yielded power (decibels) are relative to. I'm getting some arbitrary -100db-ish numbers, but I don't understand how I can have them represent a finite intensity-scale.

If I have a gradient which will represent the intensity of these numbers, where silence is black and the maximum is yellow I need to know what the maximum is. I'm really not good in the DSP-area so help would be appreciated.

For audio files you want to use dBFS, dB relative to full scale:

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


Extract feature area from image I have an image like the one above, and what I need to find are the red circled areas. They are caused by how the nerve layer is in the eye, and I have no idea what method to use to extract them. (I am looking for a brighter background from the nerve layer)
I have found all the blood vessels, so I can ignore those area easily. But I don't have any concept of how to begin the search for these areas. Note that the bright spots on the right top side (lets say 12 o clock / 1 o clock from the image center) are not caused by the nerve layer, but from being in a young persons eye.

Can someone point me to a method to extract these areas? Note that the areas could also be non-existant, or be formed like a cross, or like an hourglass.

It is going to be a tough task, and you may not achieve an automated process.

Here are some bricks that you could try:

Enhance the contrast of the image so that you can clearly distinguish all the features by looking at the image. Here is a link to an algorithm that shows good results.
In this enhanced image, extract the edges (using for example a standard Canny-Deriche edge detector).
From the edge detector output, you can start adding the required application logic in order to keep the desired vessels: size or length criteria, linear features extraction, etc.


frequency band and Bandwidth All I know so far is that frequency is how much a signal repeat it self in a given amount of time (usually one second), like shown in the following figure:

I'm confused with some concepts like frequency range or I think it's the same as frequency band (correct me if I'm wrong), so what does it exactly mean?

As I understood, if the frequency range from 1 kHz to 3 kHz, does that mean that we have a signal that maybe repeat itself  from 1000 times to 3000 time every second?

I saw that the bandwidth of analog signal is the difference between highest frequency and the lowest frequency, so in the above example is 2 kHz does that mean it is the difference between  the maximum repetition and lowest repetition of the signal?

As I understood, if the frequency range from 1 kHz to 3 kHz, does that mean that we have a signal that maybe repeat itself from 1000 times to 3000 time every second?

No.

A signal can be composed of many different repeating signals that are added up. So, if you have a sine that repeats 1000 times in a second, and a sine that repeats 3000 times in a second, and you add these two up, you get a signal with energy at 1000 Hz and 3000 Hz.

But not all signals are purely sinusoidal. Those have energy not only at discrete frequencies, but over a range of frequencies.

There's different definitions of it, but essentially, the bandwidth of a signal is the width of the frequency range in which the signal has energy.

The way your material looks, you'll pretty soon learn about the Fourier Transform and  what Power Spectral Density is. That explains the mathematical foundation of how signals can have a continuous distribution of energy over a range of frequencies instead of just at discrete (i.e. separate, countable) frequencies.


Poisson summation formula and periodic summation of Fourier transforms One of the forms of the Poisson summation formula is

$$
\sum_{n=-\infty}^{\infty} T\cdot x(nT)\ e^{-i 2\pi f T n}\;
{=} \;
 \sum_{k=-\infty}^{\infty} X\left(f - k/T\right),$$

where $x(nT)$ are samples of a continuous function $x(t)$, and $X(f)$ the fourier transform of $x(t)$. The RHS is a periodization of $X(f)$.

What I couldn't understand is the following, we can define two functions $x_{1}(t)$ and $x_{2}(t)$ that are equal at the sample points, $x_{1}(nT)=x_{2}(nT)$. We  have
 $$
\sum_{n=-\infty}^{\infty} T\cdot x_{1}(nT)\ e^{-i 2\pi f T n}\;
{=} \;
 \sum_{k=-\infty}^{\infty} X_{1}\left(f - k/T\right),$$
and
 $$
\sum_{n=-\infty}^{\infty} T\cdot x_{2}(nT)\ e^{-i 2\pi f T n}\;
{=} \;
 \sum_{k=-\infty}^{\infty} X_{2}\left(f - k/T\right),$$

where $X_{1}(f)$ and $X_{2}(f)$ are the Fourier transforms of $x_{1}(t)$ and $x_{2}(t)$ respectively.

But since $x_{1}(nT)=x_{2}(nT)$, we then also have $$
\sum_{n=-\infty}^{\infty} T\cdot x_{1}(nT)\ e^{-i 2\pi f T n}=\sum_{n=-\infty}^{\infty} T\cdot x_{2}(nT)\ e^{-i 2\pi f T n},$$

Which means $$\boxed{\sum_{k=-\infty}^{\infty} X_{1}\left(f - k/T\right){=} \;
 \sum_{k=-\infty}^{\infty} X_{2}\left(f - k/T\right)}$$

But we can choose the functions $x_{1}(t)$ and $x_{2}(t)$, even though they are equal at the sample points, in such a way that their transforms $X_{1}(f)$ and $X_{2}(f)$ are quite different from each other, and, therefore, so will their periodizations. And so why should this last equation be true in general?

But we can choose the functions $x_1(t)$ and $x_2(t)$, even though they are equal at the sample points, in such a way that their transforms $X_1(f)$ and $X_2(f)$ are quite different from each other, and, therefore, so will their periodizations.

That's a wrong conclusion. If $x_1(nT)=x_2(nT)$ $\forall n$ then the discrete-time Fourier transforms of the respective sequences must be identical, and so must be their periodized spectra.

Note that $x_1(nT)=x_2(nT)$ and $x_1(t)\neq x_2(t)$ implies that at least one of the two functions is not sampled according to the sampling theorem, i.e., the function is not completely represented by its samples. Because otherwise equal samples would imply equal functions. This might be the source of your confusion. If both functions are band-limited and if the sampling theorem is satisfied then $x_1(nT)=x_2(nT)$ implies $x_1(t)=x_2(t)$.


How to calculate baudrate? I'm going through this article https://www.elttam.com.au/blog/intro-sdr-and-rf-analysis/ and I find it hard to understand the following part:

  However, in our case we might want to transmit this signal later and need to represent each symbol with 1’s and 0’s, for instance the short pulse symbol may be represented by 0001. For this, we will measure the shortest pulse or gap that appears. Using the shortest pulse, we get a baud of 2960.

How did he calculate this 2960 rate? And why 0 should be represented as 0001 and 1 as 0111 - why 4 bits?

Thanks for your help

The article mentions that the sequence is likely pulse duration modulation, but the RfCat library doesn't support PDM. Instead, it supports OOK. So what the author is doing is representing PDM as OOK.

In particular, the author has likely measured the length of the pulses and determined that the short pulse is 1/3 as long as the long pulse, and that the entire section is 4x as long as the short pulse. So that means that OOK with 4 symbols can be used to represent 1 PDM symbol. The OOK sequence {0, 0, 0, 1} represents the PDM symbol 0, and {0, 1, 1, 1} represents the PDM symbol 1.

The baud rate here is taken as the OOK equivalent baud rate. Measuring that is just a matter of using the baudline tool there to select the region in time where the short pulse starts and stops and measuring how long that section is. That gives you a duration. You can invert a symbol duration to get a baud rate. 1/(2960 Hz) is about 0.337 milliseconds.


Notch Filter - Bandwidth vs pole / zero location vs Q factor I am reading up about notch filters. I understand the closer the pole is to the zero, the narrower the width gets and thus the bandwidth gets smaller. Is there an equation that links the bandwidth to the location of a pole? I have tried searching for so long and I couldn't get any useful answers.

I want to be able to determine the position of the pole given the required bandwidth and sampling rate and the frequency of the noise. I do not want to do it via trial and error

Thanks

In this paper by Nehora, "A Minimal Parameter Adaptive Notch Filter With
Constrained Poles and Zeros" IEEE Trans ASSP, vol 33, No 4, 1985 link

The Bandwidth is given as
$$BW=\pi(1-\rho) $$
where $\rho$ is the radius of the pole. The zero is assumed to be on the unit circle.

Excerpt from the relevant section of that paper:

And the references 1 through 3:

Another paper gives the same bandwidth and it is listed as the 3dB bandwidth - see (Yong, Kechu, Xiaofan, "A gradient algorithm for complex adaptive notch filters" - Intl Conf on Signal Processing, 2004) link.

Yet another paper with bandwidth details is by Regalia,"A Complex Adaptive Notch Filter" IEEE Sig Proc Letters, Vol 17, No 11, 2010, Link.
Note that Regalia's notch filter is based on All-pass filters and so it has a different equation for bandwidth.

The relevant section from the Regalia paper:


Can addition of two band limited signals create aliasing? Can mixing/adding two band-limited signals create any frequencies above Nyquist?

Ideally no, adding 2 signals should result as a simple sum.

provided you don’t clip

DSP systems have finite dynamic range. A double float has a very large dynamic range so it is not usually a consideration.

Fixed point arithmetic has uniform quantization but the finite dynamic range requires attention.

Floating point is not a uniform quantization, so arithmetic also introduce rounding, which are discontinuous on small scales.


Probability distribution of windowed cross-correlation This question is in the context of time-delay estimation. Say I have a stationary Gaussian stochastic process $g$, and I know its autocorrelation function $R_g(\tau)$. To do time-delay estimation, I'm computing a windowed cross correlation between $g$ and a delayed version of it. In other words,
$$
g_1 = g(x-D) \\
\phi(\tau) = \int_{-T/2}^{T/2} g(x) g_1(x + \tau)
$$
and I'm going to determine the delay by finding the maximum of $\phi$.

My question is, is it possible to get an expression for the probability distribution of $\phi$?

Expanding on my comment, $\{\phi(\tau)\}$ is a non-stationary nonGaussian random process, and I doubt that there is any simple answer (or even a rather complicated one) for the probability density function of the random variable
$\phi(\tau)$ for an arbitrary value of $\tau$.  But, the (time-varying) mean
function of the process is easy to calculate.  We have

$$\begin{align*}
E[\phi(\tau)] &= E\left[\int_{-T/2}^{T/2} g(t)g(t-D+\tau)\, \mathrm dt\right]\\
&= \int_{-T/2}^{T/2} E[g(t)g(t-D+\tau)]\, \mathrm dt\\
&= \int_{-T/2}^{T/2} R_g(\tau-D)\, \mathrm dt\\
&= T\cdot R_g(\tau-D)
\end{align*}$$

where $R_g(\cdot)$ is the autocorrelation function of the input process
$\{g(t)\}$.  Note that it is not necessary that the input process
be Gaussian for this to hold; wide-sense-stationarity is enough.
Since autocorrelation functions have a peak at the origin,
we see that $\phi(D)$ has
the largest mean value.  Also, the mean value decays away symmetrically
about $D$:  that is, $E[\phi(D+\epsilon)] = E[\phi(D-\epsilon)]$
and
$$|E[\phi(\tau)]| \leq E[\phi(D)] = T\cdot R_g(0).$$

Finding the variance of $\phi(\tau)$ is a much messier calculation
that may or may not be included in the paper cited by
Charna (which is behind a paywall).


How to classify audio files based on their ZCR, energy, and enthropy of energy features? Let's say I wanna classify the audio files of race car exhaust and normal car exhaust. I attain few features such as ZCR, energy and enthropy of energy from  the training set and I wanna build the prediction model out of those. How do I do that? Is there some machine learning model that I could use?

This is a machine learning problem, and it can be classified using any good classification algorithm like 1) Linear Classifiers: Logistic Regression, 2) Naive Bayes Classifier, 3) Nearest Neighbor, 4) Support Vector Machines, 5) Decision Trees, 6) Boosted Trees, 7) Random Forest, 8) Neural Networks. Matlab has in-built classifiers, and you can train your model very easily.

As you have already found out the features you should get some good classification results using any of these classifiers. All of the features that you mentioned are time domain-based. Also, you can incorporate some frequency-based features like 1) Spectral Centroid and Spread (C.G and deviation from CG), 2) Spectral Entropy, 3) Spectral Flux (Measure of Spectral Change over successive frames), 4) Spectral Roll-off (Spectral shape descriptor of an audio signal) which may increase the classification accuracy.

For more audio signal features you can go through this book "Giannakopoulos, Theodoros, and Aggelos Pikrakis. Introduction to audio analysis: a MATLAB® approach. Academic Press, 2014."


Digital filter selection for sensor I am having a LM35 sensor which gives $10\textrm{ mV}$ signal per degree rise of temperature. Now I have one question, which digital filter will be the best for this ?

This sensor will be use in very noisy industrial environment. Also digital filter have some cutoff frequency what should be the cutoff frequency set for this filter ? As signal coming from sensor is an analog signal but not a periodic frequency signal.

One I know is FIR Moving average filer. This takes the mean of previous samples. This I do not want to use as I have used it before,
$$
\frac 1N\sum_{n=0}^{N-1}x[k-n]
$$

If I am right then DFT is not required as I am not interested in the spectrum analysis of the signal coming from LM35 sensor.

Is there any other DSP digital filter which I can use ?

It depends a lot on the kind of noise and how fast you need the sensor to track the true temperature, and how fast you need to read the temp.

If you expect large outliers, a median filter is robust to those kinds of perturbations and they can have a relatively short effective window, and a corresponding fast response.

Following the median filter with a conventional low pass is used in some applications.

As was pointed out in the comments, the better you can define your requirements, the better you can make your measurements.

The nice thing about using a median,it is a rank statistic that it is core to nonparametric statistics, where fewer assumptions are made by the underlying probability distribution and there is a large literature on how to design hypothesis tests.


Selecting a Digital Signal Processor for a Specific Image Sensor I am working on developing a camera and I am trying to select hardware for the camera. I believe that I may need a digital signal processor to compress image data before it is transmitted to a computer or other device (MCU) . What I am trying to figure out is how to estimate the processing power needed given a specific DSP. For example I was considering the OV5670 Image Sensor. This sensor is capable of transmitting 5 megapixel image data at 30 fps. If each pixel is stored as a 10 bit value that means in total 1.5Gbps would be being transmitted to a digital signal processor of my choice. The reason I am doing compression is that my ideally the camera would record data and then upload it to a cloud based server. I think it would be ideal then for compression to take place before uploading the data to reduce bandwidth usage.
From here I am uncertain how to proceed. For example I examined the ADSP-BF529BCPZ DSP chip. This chip can execute instructions at a rate of 400 MHz. In the data specs it is specified that the processing core of this chip contains
"Blackfin processor core contains two 16-bit multipliers, two 40-bit accumulators, two 40-bit ALUs, four video ALUs, and a 40-bit shifter."
Now each MAC can do a 16 bit multiplication and accumulate the result in the 40 bit shifter. There are 2 such MACs. If I were trying to compress the data by taking a DFT and then selecting the largest coefficients then I could say to calculate a single coefficient in the frequency domain I would need to multiply each pixel by a value and then sum it. Between the 2 MACs this would take 2.5 million cycles for one coefficient of a single frame.
Based on the above calculations it seems unlikely that the selected chip would be able to handle the compression as I need it to. But at the same time I know algorithms exists such as the FFT which can drastically boost the speed of the DFT.
What I am looking for is either confirmation that my current intuition into selecting a DSP chip is correct or alternatively that I am over/under-estimating the capabilities of DSP chips. In the later case, if possible please suggest other methods for predicting the required processing power from a DSP for a specific image compression application.
Thank You in advance, if any additional information is needed please let me know. Also feel free to suggest other solutions to my problem. For example, if you believe that I do not need a DSP and could instead achieve my goals through either a MCU or an FPGA then feel free to make such a suggestion.
edit: Also if anyone looks at the datasheets for the above listed components you will see that the image sensor transfers data via MIPI-CSI2 protocol. Which the DSP chip does not support. I am aware of that shortcoming but it is not the main driver of my current question.

Good approach to first do a rough calculation of the bandwidths you need.
Couple of remarks on that:

You forgot a factor of 3, that camera has three "color pixels" per image pixel
Your use case screams "I should be using a commercial off-the-shelf USB camera"; don't engineer something very complex if it doesn't have a value proposition over what you can buy. If the value proposition is "I can learn how to work with cameras, and DSP", I'd frankly start smaller, and learn how to work with camera images on a PC in software before moving out part after part to an embedded processing platform.
You try to downplay the camera interfacing, but it's really a critical part here: in need to interface digital buses for which their processors have no dedicated interface, engineers turn to FPGAs. You'll need one, just as the "glue" between your camera
Your approach of using an MCU betrays that you're not overly familiar with what these can and can't do. Wrong order of magnitude: for this kind of problem you need what is usually sold as "application processor", with plenty RAM, and things like interfaces for network hardware and the opportunity to run an OS reasonably fast.
You picked an oldish low-power DSP core. That 400 MHz thing can do way less than modern application processor CPUs can do. For example, every SSE-enabled x86 CPU can do two multiplications at once, just like that DSP core, but they tend to run much faster.
You're planning to stream something between 1080p and 4k video .. to the internet. You really need a fast internet uplink to do that. You might want to download a 4k video snippet and look at the amounts of bits per second that needs.

So, all in all:

You need an application processor, which take your MCU completely out of business (you will probably still need an MCU, to do things like turning on power to the different power rails of your other components in the right order, but not at all involved in handling the camera data)
You need a device with a MIPI interface. If you're able to buy in big numbers, you'd typically get a camera or smart phone System-on-Chip that has a CPU (see above) which has such an interface. In your case, it's going to be an FPGA, or a MIPI-to-USB3 bridge (cypress sells these). In any way you do it, getting the interface right isn't trivial.
I don't really see how a DSP core fits in here in practice: for the data / camera handling, your application CPU does the most. For the video encoding, you'd typically use a dedicated encoding accelerator that does the computationally hardest part of it.

All in all, this sounds like you'd start with a small x86 mainboard (there's way smaller form factors than µATX, but honestly, µATX isn't so large for prototyping), and get a camera that you can attach to that - and that would usually simply be a high-quality webcam or DSLR.
For highly integrated systems, going the FPGA route is very attractive, and you can do a lot of computation by writing the right hardware description - problem is, that takes a lot of experience. For even higher integration, you'll find ASICs that are both pretty camera-specific and use-case-specific that take on the task of handling everything from camera to radio interface. In fact, DJI has bought a semiconductor manufacturer just so that they have that kind knowledge in-house, although they are by no means a semiconductor company – they sell video drones.
Honestly, I'd call the FPGA approach "unrealistic within the next year", given your current understanding of the matter, sorry. The ASIC route is impossible, since production of a single high-speed wafer costs millions, and you don't have a silicon design team.


Statistic that is responsive to changes in time series, yet not too volatile Let's say I have a fairly volatile time series $X_t$ - it doesn't have any reason to show an upward / downward trend, but it does show drops and spikes from time to time. It can also change level (e.g. for a few weeks it will be at one level, but it will quickly drop down to a new level a few % points down and stay there).
I'd like to have a statistic (basically a time series Y_t) that signals a change in $X_t$ quickly. It should also converge to the same (or close to the same) level as $X_t$ with as little delay as possible. Finally, $Y_t$ should have as little volatility as possible (I'm aware there's a tradeoff between low volatility and responsiveness).
I've tried the rolling mean - if there's a short duration (say 1 period) drop in $X_t$, then the rolling mean reacts quickly to the drop, but shows a delay in reacting to the recovery in $X_t$ (i.e. when $X_t$ is restored to normal values). The delay is the same as the window size minus the drop duration.
At this stage, I'm considering two statistics. The exponential moving average with a window of $5$:
$$E_t=\frac{\sum_{i=0}^t(1-\alpha)^iX_{t-i}}{\sum_{i=0}^t(1-\alpha)^i}$$
where $$\alpha=\frac{2}{\text{window}+1}$$
The other statistic I'm considering is a combination of EMA and median. Let $M_t$ denote the rolling median with a window of $9$. Then
$$C_t=aE_t+(1-a)M_t$$
To test their response to a drop, I did a quick simulation as follows:
a = pd.Series(np.hstack((np.ones(10)*0.81, np.ones(4)*0.77,
              np.ones(4)*0.72, np.ones(5)*0.77, np.ones(10)*0.81)))
x = 0.5
b = a.ewm(span=5).mean()
c = x*a.ewm(span=4).mean() + (1-x)*a.rolling(7).median()

figure(num=None, figsize=(27, 9))
plt.plot(a, marker="o")
plt.plot(b, marker=".")
plt.plot(c, marker=".")

By the time $C_t$ converges to the same values as $X_t$, the difference between $C_t$ and $E_t$ is really small (less than a percentage point).
$C_t$ has the feature that it won't drop as sharply as $E_t$ in response to a drop in $X_t$, but after a few periods it will abruptly drop down to the level of $X_t$. Similarly when $X_t$ recovers, $C_t$ initially shows a slower rise than $E_t$, but then suddenly shoots up, overtaking $E_t$ and catching up with $X_t$.
Now intuitively it seems that overall the EMA $E_t$ would be a better statistic, but I'm not sure. Is there any concrete, rigorous justification for the superiority of any one statistic over the other? For example, is it possible to show that maybe $C_t$ reacts similarly to two different scenarios, and is therefore unsuitable to differentiate between those scenarios?
Finally, given the requirements of quick response, quick convergence and low volatility, is there any other statistic I may have missed that's better than these two?

Hi: I don't think that there are general statements that can be made because the median component used in the second approach is pretty ad-hoc.

You probably already know this but the $\alpha$ that you defined for the EWMA is such that your EWMA has the same weighted age as just taking a moving average of length "window". Therefore, if you an moving average of length window and an EWMA with parameter $alpha = \frac{2}{window + 1)$, they will look quite similar.

The way you look at things in terms of "reaction time" till convergence or catching up is interesting.

The standard approach for this ( and more systematized but I'm not critcizing your creative approach ) is to just use an ewma with parameter $\alpha$.
The different values of $\alpha$ chosen will equate to the half life of an impulse response for $S_t$. So, just by plotting the impulse response, you can see what $\alpha$ you prefer.

So, for example, suppose you used the EWMA with parameter $\alpha$. ( $\alpha$ is between zero and one for invertibility purposes ) Therefore,
your exponentially smoothed value, $\hat{S}_{t}$, at every time step, $t$ would be:

$\hat{S}_t = \alpha \times \hat{S}_{t-1} + (1-\alpha) \times X_{t} $.

Now,  we assume that $S_{0} =1 $ and $X_{t} = 0 $ for $t$ greater than zero and that we start start smoothing $S$ at $t=1$. So, since  the halflife corresponds to the period when $\hat{S}_t$ becomes half of what it was originally, this corresponds to

$\hat{S}_{halflife} = \frac{1}{2}$.

Therefore, we have $\frac{1}{2} = \alpha^{halflife} $.

Taking the log of both sides gives $halflife = \frac{log(1/2)}{log(\alpha)}$.

So, the point is that, given the impulse response at time $t=0$, then different $\alpha$ correspond to different halflives
for the response and you can play around with the halflife to see what works best for your problem.

Finally, if this was obvious or you were already familiar with it, my apologies. Also, IMHO, the best textbook exposition of this is Brown's "Smoothing, Forecasting and Prediction of Discrete Time Series". There are used versions on Amazon that are pretty cheap.


Does windowing affect Parseval's theorem? I have computed Fourier transform from ecg data.It does obey Parseval's theorem,
relation

$$\sum_{n=0}^{N-1} \Big| x[n] \Big|^2  = \frac{1}{N}  \sum_{k=0}^{N-1} \Big| X[k] \Big|^2$$

is fulfiled. But after I use window function (Blackman,Hann,Hanning) energy is not conserved anymore.Is this normal? If not,how can I normalize my windowed fft?

I calculate fft this way:

    N = signal.shape[0]
    frq = np.fft.fftfreq(N,freq)
    frq = frq[range(N/2)]
    window = window_func(N)
    fft_y = np.fft.fft(signal*window)/np.sqrt(N)
    fft_y = Y[range(N/2)]
    plt.plot(frq, abs(fft_y))

And thats how it look when I plot. I'm a little worried about how almost all windows drastically reduce amplitude,again is this normal?

Thank you in advance for all help.

Parseval's theorem will hold, but take into account that your signal in the time domain will no longer be $x[n]$. Namely, if you have that

$$\sum_{n=0}^{N-1} \Big| x[n] \Big|^2  = \frac{1}{N}  \sum_{k=0}^{N-1} \Big| X[k] \Big|^2$$

then, if you window the signal $x[n]$ with a window $w[n]$, your signal will now be $\hat{x}[n]=x[n]w[n]$, and the theorem will now return the following equality:

$$\sum_{n=0}^{N-1} \Big| \hat{x}[n] \Big|^2  = \frac{1}{N}  \sum_{k=0}^{N-1} \Big| \hat{X}[k] \Big|^2$$

where $\hat{X}[k]$ is the DFT of $\hat{x}[n]$ (and obviously differs from $X[k]$).


How to design very narrow, sharp lowpass filters - Only DC needed? I asked exact the same here one year ago: One year ago question I asked this question again here, because after one year, I have learned more about DSP and am thinking ways to improve it.

Filter description

Filter input data rate is 64kHz
Filter output data rate is 0.5Hz - 2Hz
The filter is preferably to have bandwidth of 1mHz with very sharp transition area, since only DC is needed and the noises at 0-5Hz is large.

The answer I got

The averaging data is the best DC I can get. I bought this answer. However, it now becomes plausible to me, since a filter always have a bandwidth and frequency responses. Averaging is essentially a FIR with no superior characteristics.

Questions

Shall the filter be designed in this way: 9 stages of halfband filters(64kHz -> 125Hz). Followed by three stages of cascaded integrator–comb (CIC) filters with R=5. Then the output will be fed into a 1st order IIR filter?
Maybe 16 stages of halfband filters? (64kHz ->  0.9765625Hz) And then be fed into an IIR filter.
How shall the filter be designed?

Your questions still leave me wondering as to what you're actually designing. For software implementation on modern x86 CPUs, CICs make almost no sense, but they are extremely elegant in hardware.

These filter definitions are ridiculous if you're planning to use a FIR – a transition width of 1mHz means that a minimum phase equiripple filter [1, (5.75), p. 142] would have

$$
N_e' \approx \frac 16 \log_{10}\left(\frac{1}{10\delta_1\delta_2^2}\right)\frac{f_s}{\Delta f}
$$
taps, with $f_s$ the sampling rate, $\Delta f$ the transition width, and $\delta_{1,2}$ the passband ripple and the stopband attenuation, respectively.

Now, your use case; I'm guessing from your the graphs from your question from 1 year ago, and also bearing in mind that you'd need incredible numerical accuracy to make use of > 80dB attenuation:

$$
\begin{align*}
f_s &= 6.4\cdot10^4\,\text{Hz,}\\
\Delta f &= 10^{-3}\,\text{Hz,}\\
\delta_2 &= -90\,\text{dB}\\
&= 10^{-9}\\
\delta_1 &= \delta_1\,\text{, everything else being wasteful either on suppression or precision side,}\\\\
N_e' & \approx \frac 16 \log_{10}\left(\frac{1}{10\cdot10^{-27}}\right)\frac{6.4\cdot10^4}{10^{-3}}\\
&=\frac16\cdot 26\cdot6.4\cdot10^7\\
&\approx 2.7 \cdot 10^8\,\text.
\end{align*}
$$

Obviously, Bellanger's approximation wasn't meant for such extreme cases, but the order of magnitude alone highlights the impossibility. Of course, you'd implement this as a decimating polyphase FIR, which means that you'd only need to apply
$$\frac{N_e'}{\frac{f_{s,input}}{f_{s,output}}}=N_e' \frac{f_{s,output}}{f_{s,input}}= N_e'\frac{\frac12}{6.4\cdot10^4}\approx2110$$
tapped filter every input sample clock, but the filter length puts impossible constraints on your machine accuracy if you're using floating point, and very hard constraints if you're building this in hardware on the bitwidth of your samples. As you notice, it does make a lot of difference for what we're designing a filter here – Matlab is not hardware :)

Yes, you can split your filter into multiple stages, but that will usually only reduce the computational load by running "parts" of the filter at lower rates, and will not solve the precision/bit width problems.

So we're ditching FIRs here, unless you're willing to reduce the requirements on stopband attenuation or transition width.

On to IIRs:

There's an approximate formula for the order necessary to achieve given specifications, and that is [1, (7.18), p. 184]:

$$
N_{IIR}\approx 1.076\log_{10}\left(\frac{2}{\sqrt{\delta_1}\delta_2}\right)
\log_{10}\left(\frac{f_s}{\Delta f}\frac4\pi\sin\left(2\pi\frac{f_1}{f_s}\right)\right)\text,
$$
$f_1=f_2-\Delta f$ and $f_2$ being your pass- and stopband edge. Since you're lowest output rate is 0.5 Hz, I'm guessing $f_2=0.5\,\text{Hz}$.

Applying the same specifications as above,
$$\begin{align*}
N_{IIR}&\approx 1.076\log_{10}\left(\frac{2}{\sqrt{10^{-9}}10^{-9}}\right)
\log_{10}\left(\frac{6.4\cdot10^4}{10^-3}\frac4\pi\sin\left(2\pi\frac{4.99\cdot10^{-1}}{6.4\cdot 10^{4}}\right)\right)\\
&\approx 1.076\log_{10}\left(2\cdot10^{13.5}\right)
\log_{10}\left(\frac{2.56\cdot10^8}\pi\cdot4.899\cdot10^{-5}\right)\\
&\approx 1.076\log_{10}\left(2\cdot10^{13.5}\right)
\log_{10}\left(3992\right)\\
&\approx 1.076\cdot13.80\cdot3.601\\
&\approx 53.48
\end{align*}
$$

Note that this is but an estimate! Also notice that this is far from a design approach, it just gives us an idea of how complicated the IIR will end up being.

Now, realistically, let's round up to $N_{IIR}=64$, and then do the classical approach of directly calculation:

We'll consider elliptic filters, which are optimal in the sense that for any given order, and accepting a certain ripple (which is your DC estimate accuracy, here), they have the most narrow transition width [1, p. 183].

Now, designing these taps isn't a trivial feat; scipy gives me something that roughly looks like that:

Notice the enormous dynamic range of each tap set ($\approx 10^{18}$), and of the overall filter ($\approx10^{65}$). That filter is nothing that reaches its accuracy on a double precision IEEE754 float.

I propose a hybrid approach: take cascaded halfband filters to get down to something like $\frac1{512}$ of the input sample rate, and then apply an elliptic filter.

1 Bellanger, Maurice: Digital Processing of Signals: theory and practice, second edition, Wiley 1988.


From pressure to SPL : microphone calibration using sound level meter calibrator and sound level meter I have a microphone + signal conditioner with adjustable gain and filter. I capture amplified signal using oscilloscope. I'm using this signal to compute equivalent continuous sound level

With signal conditioner, I set needed gain and A-weighted filter.
From oscilloscope signal V I can get pressure p = V/ (gain x mic sensitivity)

We can suppose A-weighted filter filters out frequency out of audio range. So I can easily compute leq = 20 log10(rms(pa)/pref), pa is time domain pressure A filtered, pref is reference pressure = 20e-6 pa

I used a sound level calibrator to check if it works. I get a leq of ~94dB and 114dB when applying 94dB and 114dB level with calibrator.

I also have a sound level meter, I also checked it with sound level calibrator.
Now, in my room (some noise because of oscilloscope...) with A filter it measures ~46 dB(A). Now I put my microphone close to it, with microphone I get a leq of 55 dB. When changing sound level meter filter from A to C results are better : sound level meter measures ~56dB.
I changed microphone acquisition duration but leq results are same.

How can we explain such leq difference when sound level meter and
microphone measure dB(A)?
How can I calibrate my microphone with    sound level meter?

First, I would check the positioning of sound level meter and microphone. If you are putting them close to the oscilloscope a few centimeters of displacement may change the measured value dramatically. (6dB/doubling of distance)

Second, which microphone do you use? The link you provided leads to a preamplifier. The actual capsule can either be a free field or a pressure field one. A free field capsule should be facing the sound source directly, while a pressure field capsule should face it sideways.


Determining Loss of Information by Taking Average (Mean of Signal) So basically information is defined by expected value of Shannon's information i.e. Entropy. I am curious how much information is lost if we simply take the average of the sample given to us. I am also curious if we also take the second moment (variance) of the same signal, how much information do we gain by using the variance along with the mean? In other words:

$$H_{\textrm{Mean}}(\Sigma_0^n x_iP_i)=\ ?$$
$$H_{\textrm{Variance}}(\Sigma_0^n (x_i-\mu )^2)=\ ?$$

Motivation: I just want to know how much information is dropped if we simply drop variace term and use just mean instead. I believe this concept cant be expanded to also consider information in third and fourth moments.

Consider observing $y(n) = w(n)$ where $w(n)$ are IID gausian distributed with mean 0 and variance $\sigma^2$. The entropy of each other observations $y(n)$ is thus: $$\log_22\pi e\sigma^2$$.

Now if you average $y(n)$ the resultant variable is also gaussian distributed with mean 0 but variance $\sigma^2/N$. Therefore the entropy is $$\log_22\pi e\sigma^2/N$$ so the entropy difference is $$\log_2N$$ and this would be the entropy reduction.

The amount of information lost will depend on the distribution of the random variable. However, there are some interesting problems in signal processing where there would essentially be no loss.

For ex: consider a set of observations taken as $$y(n) = x + w(n)$$ where $w(n)$ is normally distributed IID random varibles and suppose we take $N$ observations of this kind.

Then the mean of $y(n)$ essentially represents what we refer to as a sufficient statistic for $x$ and infact represents the minimum variance unbiased estimator, i.e. the most efficient estimator. So in this case the sample mean encompasses all the necessary information and no information is lost.

However, if the noise $w(n)$ were  not gaussian we couldn't claim this. So as you see the amount of information lost by sample mean depends on the distribution of the underlying variable.

Another example would be a trivial one of a deterministic random varaible. Here also no information would be lost by the sample mean.

Another way to look at it would be if the likelihood of posterior distribution (parameterized by the sample mean), after observing the sample mean renders the posterior useless, meaning encompasses no additional information about the random variable then there is no information loss, otherwise we would incur information loss.

Hope that helps


Discrete-time Convolution Convergence Issue I'm trying to solve a problem on convolution from Alan V.Oppenheim:
Find the convolution output $y[n]$ for the following signals:
$$x[n]= u[n]\quad\text{and}\quad h[n]=a^{-n}u[-n-1], \ a>1 $$
I started the evaluation:
$$y[n]=\sum_{k=-\infty}^{+\infty} u[k]a^{-n+k}u[-n+k-1]$$
considering that $u[k]=1$ for $k\ge0$ and $u[-n+k-1]=1$ for $k\ge n+1$
which I evaluated to
$$y[n]=a^{-n}\sum_{k=n+1}^{+\infty} a^k$$ where $n$ could be $<0$ or $>0$
and I tried to evaluate for $n>0$:
where I faced an issue:
I subsittuted $k-1$ to $m$ and reframed the equation as:
$$y[n]=a^{-n+1}\sum_{m=n}^{+\infty} a^m=a^{-n+1}\left[\left(\sum_{m=0}^{+\infty} a^m\right)-\left(\sum_{m=0}^{n-1} a^m\right)\right]$$
now for $a>1, \ \displaystyle \sum_{m=0}^{+\infty} a^m$ will not converge. How will I evaluate this?

In the third edition of Oppenheim and Schafer's Discrete-Time Signal Processing, the corresponding exercise is given with $h[n]=a^nu[-n-1]$, $a>1$ (note the positive sign in the power of $a$). In that case everything works out fine. So either there's a typo in your edition or you made a mistake copying the exercise.


What Is the Walsh Hadamard Transform? How Could One Use It in Image Processing? I am trying to teach myself about the WHT but there dont seem to be many good explanations of it online anywhere. I think I have figured out how to calculate the WHT, but I am really trying to understand why it is considered useful within the image recognition domain.

What is so special about it, and what properties does it bring out in a signal that would not show up on classical Fourier transforms, or other wavelet transforms? Why is it useful for object recognition as pointed out here?

NASA used to use the Hadamard transform as a basis for compressing photographs from interplanetary probes during the 1960's and early '70s.  Hadamard is a computationally simpler substitute for the Fourier transform, since it requires no multiplication or division operations (all factors are plus or minus one).  Multiply and divide operations were extremely time intensive on the small computers used on board those spacecraft, so avoiding them was beneficial both in terms of compute time and energy consumption.  But since the development of faster computers incorporating single-cycle multipliers, and perfection of newer algorithms such as the Fast Fourier Transform, as well as the development of JPEG, MPEG, and other image compression, I believe Hadamard has fallen out of use.  However, I understand it may be staging a comeback for use in quantum computing.
(NASA use is from an old article in NASA Tech Briefs; exact attribution unavailable.)


Finding zebra-like pattern in image (Detection of structured-light fringe centerline from photo) I'm working in a project where fringes are projected against a subject, and a photo is taken. The task is to find the centerlines of the fringes, which represent, mathematically, the 3D curve of intersection between the fringe plane and the subject surface.

The photo is a PNG (RGB), and former attempts used grayscaling then difference thresholding to get a black-and-white, "zebra-like" photography, from which it was easy to find the midpoint of each pixel column of each fringe. The problem is that, by thresholding and also by taking the mean height of a discrete pixel column, we're having some precision loss and quantization, which is not desired at all.

My impression,  by looking at the images, is that the centerlines could be more continuous (more points) and smoother (not quantized) if they were detected directly from the non-thresholded image (either RGB or grayscale), by some statistical sweeping method (some flooding / iterative convolution, whatever).

Below is an actual sample image:

Any suggestion would be much appreciated!

I suggest the following steps:

Find a threshold to separate the foreground from the background.
For each blob in the binary image (one zebra stripe), for each x, find the weighted center (by pixel intensity) in y direction.
Possibly, smooth the y values, to remove noise.
Connect the (x,y) points by fitting some kind of curve. This article might help you. You can also fit a high-level polynomial, though it is worse, in my opinion.

Here is a Matlab code that shows steps 1,2 and 4. I skipped the automatic threshold selection. Instead I chose manual th=40:

These are the curves that are found by finding the weighted average per column:

These are the curves after fitting a polynomial:

Here is the code:

function Zebra()
    im = imread('http://i.stack.imgur.com/m0sy7.png');
    im = uint8(mean(im,3));

    th = 40;
    imBinary = im>th;
    imBinary = imclose(imBinary,strel('disk',2));
    % figure;imshow(imBinary);
    labels = logical(imBinary);
    props =regionprops(labels,im,'Image','Area','BoundingBox');

    figure(1);imshow(im .* uint8(imBinary));
    figure(2);imshow(im .* uint8(imBinary));

    for i=1:numel(props)
        %Ignore small ones
        if props(i).Area < 10
            continue
        end
        %Find weighted centroids
        boundingBox = props(i).BoundingBox;
        ul = boundingBox(1:2)+0.5;
        wh = boundingBox(3:4);
        clipped = im( ul(2): (ul(2)+wh(2)-1), ul(1): (ul(1)+wh(1)-1) );
        imClip = double(props(i).Image) .* double(clipped);
        rows = transpose( 1:size(imClip,1) );
        %Weighted calculation
        weightedRows  = sum(bsxfun(@times, imClip, rows),1) ./ sum(imClip,1);
        %Calculate x,y
        x = ( 1:numel(weightedRows) ) + ul(1) - 1;
        y = ( weightedRows ) + ul(2) - 1;
        figure(1);
        hold on;plot(x,y,'b','LineWidth',2);
        try %#ok<TRYNC>
            figure(2);
            [xo,yo] = FitCurveByPolynom(x,y);
            hold on;plot(xo,yo,'g','LineWidth',2);
        end
        linkaxes( cell2mat(get(get(0,'Children'),'Children')) )
    end
end

function [xo,yo] = FitCurveByPolynom(x,y)
   p = polyfit(x,y,15);
   yo = polyval(p,x);
   xo = x;
end


Measuring frequency response with white noise signal I need to measure a frequency response of an eleсtric circuit (simple like an RC-filter). I suppose that I can generate an input signal that has spectrum like white noise (for example, m-sequence signal or sine sweep), measure an output signal of a circuit and do calculations. And I can get a spectrum if I divide an output signal spectrum by an input signal spectrum.
Am I right? Or am I describing a wrong method?

The question is about measurements in real world, not about simulations, am I right? If this is the case, there is an established technique and instrumentation for doing such measurements. The generic name of an instrument commonly used in laboratories and field testing for this purpose is Vector Network Analyzer. The article describes the operation principle in details and contains references for further research into the subject.
There is another kind of applications (audio environment acoustics, power grid control) where you need to measure frequency response of the system in real time (on the fly) -- or, admittedly, to do relatively instantaneous measurements of rapidly changing frequency response. In measurement of this kind, if the device under test (DUT) is linear, you can leverage the impulse response/frequency response duality, implementing pulse testing of your DUT. To be able to test a non-linear DUT, you can lower the crest factor (peak power/average power index) using a pseudorandom noise (as maximum length sequences that you mention) instead of single pulse stimulus.
But, to measure frequency response of "simple" networks "like an RC-filter", you need not go into such intricacies and can do with an ordinary lab equipment.


How to narrow the frequency band of a wavelet I have a ricker wavelet with a dominant frequency of 15 Hz. The fourier transform shows its frequency band is almost to 50 Hz. How can I narrow the frequency band of this 15 Hz ricker wavelet?

I have added a picture of the wavelet and its fftshifted frequency band (The x axis is not the frequency samples).

Center frequency and bandwidth are linked in a single "scale" parameter for wavelets. That means if you change the bandwidth you also change the center frequency and vice versa.

For the discrete wavelet transform on a dyadic lattice the relative bandwidth is also fixed at approximately one octave.

If you want a continuous wavelet transform frame with a smaller relative bandwidth you need to construct a different wavelet. Intuitively, it will come with more oscillatory cycles, i.e. more zero crossings.

The Ricker wavelet you have coincides with the 2nd so called Hermite function. Higher order Hermite functions also have an increasing number of zero crossings. Unfortunately they're one of the examples of where the intuition of adding more cycles will not work. The Hermite functions grow in both frequency as well as time support with increasing order. So that's not the way to go. I'm only making this clear because you might find a reference to them in relation to the Ricker Wavelet.

So instead I would try to adapt the construction of the wavelet to your specific needs, which you will have to outline in your question to receive a more detailed answer.


Definition of "Power" in "Power Spectrum" (for translator) I'm a technical translator who recently undertook translation of documents related to signal processing. Being a novice in the field, I'm searching the web for terminology and concepts.
It appears that local web content includes two instances of translation for "power spectrum":

(electric power) spectrum
(power) spectrum

Upon my query, the requestor says it's not appropriate to translate "power spectrum" as "(electric power) spectrum". But there are quite a number of translation as "(electric power) spectrum", so I'd like to know the exact definition of "power" in the "power spectrum" to make sure that the translation as "(electric power) spectrum" is wrong.
The following is a Wikipedia page that my requestor provided as a reference:
https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density

... For continuous signals over all time, such as stationary
processes, one must rather define the power spectral density (PSD);
this describes how power of a signal or time series is distributed
over frequency, as in the simple example given previously. Here, power
can be the actual physical power, or more often, for convenience with
abstract signals, is simply identified with the squared value of the
signal.

My questions are as follows:

Is the highlighted part applicable as definition of "power" in "power spectrum" as well?

If so, is "(electric power) spectrum" incorrect translation because the "power" spans other types of physical power in addition to "electric power" (and also because "power" can imply the squared value)?
Otherwise, could someone provide better explanation?

Thank you.

the requestor says it's not appropriate to translate "power spectrum" as "(electric power) spectrum".

And I'd agree with him.
Let's start with the last of your questions:

Otherwise, could someone provide better explanation?

What's the layman area where one meets the term "spectrum"? Right, it's this rainbow:
Zátonyi Sándor, (ifj.) Fizped / CC BY-SA (https://creativecommons.org/licenses/by-sa/3.0)
What's depicted here is how a prism splits white light, which contains a continuum of frequencies of light, into components that all have a different frequency.
The power spectrum is just the answer to the question "how much power is there in the white light at the frequency of this specific color?".
No electricity involved at all. Just light.
Generally, in DSP, if you try too hard to identify magnitudes with any particular physical entitity, you'll often "phrase yourself into a corner", because DSP is just applied math, and that math doesn't care whether the signal analyzed is the amplitude of an electromagnetic wave like light, the intensity of pixels in an image, the instantaneous position of a vibrating string of a guitar, air pressure variations caused by that string, blood pressure, the magnitude of a gravity disturbance caused by rotating supermassive star pairs, or the voltage at some measurement device.
In fact, the power spectral density is as much a mathematical concept with no relation to any specific physical entity as you'll encounter: It's literally just "if you have a random function of a single real-valued number, and you can calculate the autocorrelation of that over a difference on that number axis, then the Fourier transform of that is called the power spectral density".
Physics is just good to us in that this really describes something that can be measured as power, if the random function can be understood as an amplitude. But again, nothing here says that anything is about electricity.

Is the highlighted part applicable as definition of "power" in "power spectrum" as well?

Yes. It's what I'd recommend to use. Power is just the square of some amplitude signal.

If so, is "(electric power) spectrum" incorrect translation because the "power" spans other types of physical power in addition to "electric power" (and also because "power" can imply the squared value)?

yes. and yes.


Complex impulse response functions? I think I must have misunderstood something regarding the relation between impulse and frequency responses.

I have calculated the Impulse Response Function, IRF, from the Frequency Response Function, FRF, by FFT and then I find some strange things, among these that my assumption that the IRF will always be real is not correct. I have in my code used only the real part of the result from FFT when defining the IRF. Then the amplitude response (and phase) are not correct when I convolve my signal with a sine wave. I expect no phase lag and a slight amplification of a sine wave with period 10 seconds, but get a signal with the wrong amplitude and a slight phase mismatch instead.

What should I do with the imaginary part of the IRF i get back from FFT?

Details for those with a lot of time:

I am using Python and NumPy to try to work out the response of a LTI system to an arbitrary response. I guess the same would apply to Matlab or any other numerical tool implementing DFT. I know a priori the frequency response function, FTR (called a Response Amplitude Operator, RAO, in my field, naval hydrodynamics). As a test case I define my FRF as follows:

# Define RAO
omegas = 2*pi / numpy.array([5, 9, 10, 11, 30])
rao = numpy.array([0, 0, 1.1, 0, 0], dtype=complex) # Fake "ship roll" RAO

The frequency response for my fake ship roll RAO is 1.1 for a wave period of 10 seconds and zero when you get more than one second in wave period away from the peak.

I am using numpy.fft to find the impulse response function, IRF. I interpolate my RAO and apply FFT:

def make_double(x):
    """
    Make a double spectrum by copying the input reversed
    Does not divide the results by two after the mirroring
    """
    return numpy.concatenate((x, x[::-1]))

def rao_to_irf(omegas, rao, tmax, dt):
    """
    Take an RAO defined for a set of (increasing) frequencies and
    compute the corresponding impulse response function by FFT.

    The RAO is inter-/extrapolated to the necessary delta omega
    required to obtain the wanted time step, delta t

    The RAO is assumed to be a complex array
    """
    assert omegas[0] < omegas[-1], "The RAO frequencies, omegas, must be given in increasing order"

    N = numpy.ceil(tmax/dt)
    dw = pi/N/dt

    omegas_long = numpy.arange(N, dtype=float)*dw
    rao_long = interpolate_complex_1d(omegas, rao, fill=0)(omegas_long)

    irf = ifft(make_double(rao_long))[:N]*N/tmax
    t = numpy.arange(N)*dt

    return t, irf.real

Now I have an IRF (?). To check the response I also implemented the inverse and see that the RAO/FRF calculated from the IRF looks similar to the input RAO (when interpolated to constant frequency spacing).

I can now combine my RAO/FRF with a wave spectrum and perform a time realization using FFT. I also make a realization of the spectrum to get the unmodified input wave signal. This signal I convolve with the IRF.

I get the expected result if my RAO/FRF is 1 for all frequencies (meaning I get the wave signal back). However the response to the peaked RAO defined above is not as I expect. The time realization from the wave spectrum multiplied by the RAO^2 is not the same as the result of the convolution. Should not these match. I am pretty sure our code to make realizations of response spectra works as it should and I have tested with wave spectra consisting of one pure sine wave and the response is as expected with an FFT of the response spectrum, but not when using convolution

Sorry for the convoluted question

Filters are usually real, not complex, but occasionally they are complex.  What do you do with the complex impulse response?  Convolve it with your data, just like you would with a real impulse response.  This will likely result in a complex output.  What you do with the imaginary part at that point depends on the situation.  Sometimes it makes sense to leave the data complex, sometimes it makes sense to get rid of the imaginary part.

The imaginary part of the filter would particularly affect the output phase, so the fact that the output phase is not what you expected when you only used the real part of the impulse response makes sense.


Evaluating the spectral density of generated noise through the autocovariance I'm working on generating noise signals $X(t)$ (with $t \in \left[0,T\right]$ with step size $\delta t$) with a prescribed power spectral density $S_{XX}(f)$ and I'm figuring out how well I am generating said noise; I need to know if it has the desired spectral density (in case someone notices, I posted about similar issues a while back. I've since then worked on a different topic and I am now returning to this).

However, estimating the spectral density of a signal is a tedious process. It involves a lot of estimation through periodograms and window functions and such, and it is therefore not so easy to quantify numerically if the spectrum not only has the desired trends but actually implements the desired parameters (in practice I will simply use a spectrum analyzer, but that is a different story).

The auto covariance function $C_{XX}(t') = \left<X(t) X(t+t') \right>$ on the other hand is very easy to calculate and subsequently fit. Moreover, the auto covariance function and the spectral density are related by a Fourier transform: $C_{XX}(t') = \int_0^\infty{S_{XX}(f)\cos{2\pi f}}$ (note that I'm using the cosine transform as this is customary in these settings, due to the evenness of the spectral density).

So my question is, how much do I learn about the spectral density by fitting the auto covariance? Say I am working with the Ornstein Uhlenbeck process, which has auto covariance $C_{XX}(t') = \frac{c \tau}{2}e^{-t'/\tau}$ and spectral density $S_{XX}(f) = \frac{2 c \tau^2}{1+(2\pi f \tau)^2}$. If I fit my numerically generated $X(t)$ to $C_{XX}$ to extract $c$ and $\tau$ and am able to do so with high accuracy, do I then know that my spectral density is also as desired?

Intuitively I would say yes, Fourier pairs contain the exact same information. The reason I am not sure if this is the case is due to the fact that I am unsure of how sampling rates and such come into this. Of course one cannot expect a spectral density to be accurate up to the GHz range if one has an $X(t)$ point every millisecond. But is this equally true for the auto covariance? In other words, if my autocovariance is accurate up to some specific lag time (which means I have a high enough sampling rate I suppose), is the spectral density then also accurate up to some specific frequency? And can I quantify this in some way?

What I mean by that last sentence is, can I tell up to (or down to) what frequencies I can trust my spectral density, given my auto covariance fit, $T$ and $\delta t$?

The autocovariance and spectral density are Fourier transform pairs. Therefore, the autocovariance function and spectral density are equivalent.

Have a look to the definitions and examples for white noise in this presentation.


Chromatography Baseline Placement I have been looking around on Google, and I have not been able to find an explanation of how to do baseline placement in chromatography. what are some algorithms for placing a baseline on a chromatogram, how do those algorithms work. (Referral to Python implementations are welcome).

In some physico-chemistry methods, baselines can be modelled by closed-form lineshapes. This is not typically the case for chromatography, and it can be even worse depending on the type of chromatography (gaz, liquid).

So you have to rely on some morphological properties that you believe valid for your type of signal. You can model the background by a given polynomial, assume derivability of the peaks, etc. And then fit, or optimize a cost function.

In a recent work, we have modeled gaussian peaks as sparse with sparse derivative. We add that peaks are typically positive. Then we suppose that the baseline is smooth enough to remain after the  use of a low-pass filter. Finally, we consider Gaussian noise, yet it works with Poisson noise as well.

The BEADS (Baseline Estimation And Denoising w/ Sparsity) method is described, and compared with two other methods in  "Chromatogram baseline estimation and denoising using sparsity (BEADS)", by Xiaoran Ning, Ivan W. Selesnick, Laurent Duval, in Chemometrics and Intelligent Laboratory Systems, December 2014, http://dx.doi.org/10.1016/j.chemolab.2014.09.014.

We have made the associated Matlab toolbox available in Matlab.  Contributed versions in C++, R, and Python are made available through the companion BEADS page.


Why is Prewitt filter High Pass Prewitt filters are popular filters in image processing for edge detection
http://en.wikipedia.org/wiki/Prewitt_operator

Can anyone give a proof on why Perwitt and other edge detecting filters are high pass?

That depends on the definition of high-pass filter. If you define a high-pass filter as a filter that has high response in the high frequencies in frequency domain, then the easiest way is to take a look at the magnitude of Fourier transform, (by definition).

Applying Fourier transform (in Matlab)

 A = fftshift(abs(fft2(padarray([-1 -1 -1; 0 0 0; 1 1 1],[10 10]))));

results in the following image:

Now the interpretation - The part in the middle is the low frequencies. It has low response. There are 2 high responses, both with zero X frequency and some high Y frequency. That is not surprising since  we took a filter that detects edges in Y direction.

Why does it make sense to define high-pass filter in this way?

Because a convolution can be thought as multiplication in Frequency domain. That is, if you have a signal $S$ and filter kernel $f$,

$F[S**f] = F[S] * F[f] $

Or, put in another way:

$ S**f  = F^{-1} [ F[S] * F[f] ] $

where convolution is denoted by $**$

There is another, more intuitive way, that does not involve Fourier transform. A response to a linear filter is strong when the underlying signal "looks" like the filter itself. Therefore, there should be a strong response to edges.


Relative performance of RLS and LMS filters It's known that the RLS filter converges faster than the LMS filter in general, but that if you're tracking time varying parameters the LMS algorithm can perform better.

My question is under what conditions does this hold? I understand that the LMS filter is like a point estimate, but the RLS uses more data - when would using less data be helpful?

Using less data is helpful when, as you said, the parameters are time varying, and in particular when they change a lot.

The key difference is that LMS is a Markov process.  It has its current state, but other than that it does not remember data from the past.  For time-varying signals this is a feature because past data will give you erroneous information about the current parameters.

The RLS algorithm uses all of the information, past and present, but that can be a problem if the past data is misleading for the current parameters.

If you are looking for a quantitative rule for when to use one or the other, I don't have one.  RLS is more computationally intensive than LMS, so if LMS is good enough then that is the safe one to go with.  RLS converges faster, but is more computationally intensive and has the time-varying weakness, so I would only use it if the parameters don't vary much and you really needed the fast convergence.


Question about time shifting property I know that time shifting property of discrete-time Fourier Transform states that
$x[n]\leftrightarrow X(e^{j\omega})\implies x[n-n_0]\leftrightarrow e^{-j\omega n_0}X(e^{j\omega})$
So given a signal $\;y[n]=u[n]-u[n-5]\;,\;$What is $\;\mathfrak{F}\{ y_{(2)}[n-1]\}\:?$
$ y_{(2)}[n]\;$ denotes the slowed-down version of $\:y[n].\;$And $\:y[n] \leftrightarrow e^{-j2\omega}\frac{\sin(5\omega/2)}{\sin(\omega/2)}$
From time expansion property $\:y_{(2)}[n]\leftrightarrow e^{-j4\omega}\frac{\sin(5\omega)}{sin(\omega)}\:$then follows the time shifting property
$\mathfrak{F}\{ y_{(2)}[n-1]\}=e^{-j6w}\frac{\sin(5\omega)}{\sin(\omega)}$
But the answer should be $\:e^{-j5\omega}\frac{\sin(5w)}{\sin(\omega)}$
What is wrong with my understanding?


How to find out the farrow coefficients if FIR coefficients are given? FIR filter coefficients are known. Then what is the matlab code or function that is used to determine the corresponding farrow structure coefficients?

There is no standard way to determine the coefficients of the Farrow structure for a given FIR filter. The Farrow structure is an implementation of a whole class of FIR filters with transfer function(s) $H_{\mu}(z)$ with a continuous control parameter $\mu$. Often this parameter $\mu$ determines a fractional delay, but it could as well be any other adjustable filter property.

The transfer function $H_{\mu}(z)$ of the Farrow structure is given by

$$H_{\mu}(z)=\sum_{k=0}^{K-1}C_k(z)\mu^k\tag{1}$$

where $C_k(z)$ are FIR transfer functions. From $(1)$, the impulse response is

$$h_{\mu}[n]=\sum_{k=0}^{K-1}c_k[n]\mu^k\tag{2}$$

So each filter coefficient (of the corresponding transversal filter structure) $h_{\mu}[n]$ is implemented as a weighted sum of coefficients $c_k[n]$. This is why there is no one-to-one mapping from given FIR filter coefficients to the coefficients of the Farrow structure.

You need to specify which filter property the parameter $\mu$ is supposed to control, and then you need to design the coefficients $c_k[n]$ of the Farrow structure accordingly.