Efficient calculation of exact frequency from a broad frequency range (1 Hz - 1000 Hz) I have some input signal sampled at 1-10 MS/s. Let's assume that it contains one strong frequency component (with frequency anywhere in the range from 1 Hz to 1000 Hz or even higher, changing over time) plus a number of higher harmonics (with quite a bit lower amplitude than the main frequency) and some noise.

The frequency may rise from zero to 1 kHz in a few seconds. The idea is to calculate the result approximately every 10-20 periods of the signal.

I need an algorithm which is as fast as possible to calculate the precise main frequency in real time.

I'm currently testing a couple of methods (like this or this) which can already give a pretty accurate result (under 0.01% when there are not too many higher harmonics) if I know an approximate frequency upfront and feed the FFT with cca. 10 periods (if the signal is around 100 Hz, I would use 100 ms slices; in case of 1 Hz, I'm ok with a lower accuracy, so taking just 1 - 2 periods is fine in that case). All these algorithms use just two or three bins from the FFT (around the highest amplitude) to calculate the exact frequency.

The problem is that in order to cover a wide frequency range without any idea what the base frequency is, I would first need to calculate FFT a lot of times. For example, first calculate all chunks of 10 ms just in case there is some 1 kHz component present. If it's not, maybe switch to 20 ms, then 40 ms ... and so on until we find the first signal length which adequately covers cca. 10 periods of the signal. That might slow down the algorithm 10-fold compared to the situation when frequency is approximately known upfront.

My question is: is there any way to first calculate 10 ms FFTs for each chunk of the signal, and if the frequency turns out to be too low, somehow combine results from two neighbouring 10 ms + 10 ms FFTs with a method that would require less computation resources than calculating FFT of the full 20 ms from scratch (given that I probably only need one or two FFT bins for a 20 ms signal anyway, and in this case I should know which ones).

Or any other idea of the fastest way to determine an approximate frequency quickly? (Some even use zero-crossing to do that which works fine as long as you are looking at nearly perfect sine curve, but that's not useful when higher frequencies and noise are present.)

Thanks for citing my blog article on frequency calculation.  Based on your emphasis on efficiency I am going to recommend a different approach which should get you equivalent or superior results with way fewer calculations.

The first thing you should do is smooth your signal heavily with exponential smoothing.  This will squash your harmonics and mitigate noise.  You can carry your results forward buffer by buffer so it is extremely efficient.  While you are doing the smoothing, you can also check for zero crossings.  The signal need not be a "nearly perfect sine curve" for zero crossings to give you a good estimate of the frequency.  After heavy smoothing, the zero crossing should correspond to your fundamental frequency.

If your signal is sufficiently sine like, you can calculate the frequency quite accurately by using the three point (simplest case) formula with the appropriate spacing that you can find here:  Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 1).  You will want to select the center point at the peak (or trough) and the other two points at the nearest zero crossings.  This formula will give you as good results as if you did interpolation at the zero crossings and calculated the frequency from there.  You can then calculate the frequency at every peak and trough, that is every half cycle, and generate a new sequence of values.  You can then apply a best polynomial fit on intervals of this sequence to give you very good calculations of the frequency and how it is varying.

A more accurate, but slightly more computationally expensive approach would be to calculate two bins of the DFT and use the formula found in my article: Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT.  For the best result the frequency should fall halfway between the bins so you want to select a frame interval of a whole number of cycles and a half.  Two and a half would be a good compromise between getting an accurate reading with a varying frequency and getting a better reading with multiple cycles.  Point spacing is key for saving calculations.  Four points per cycle will give the best results, so you want to be close to that.  For two and a half cycles, that means ten points, which is not many.  Since you are calculating the bin values separately, there is no advantage to a power of two.  You can build lookup tables for the ten points for the two bin and three bin ahead of time to make it even more efficient.  Pick a uniform spacing parameter that gets you as close to two and a half cycles as you can.  Your calculations result in a cycles per frame value, then use your spacing value to convert to cycles per whatever your measure is.

Feel free to add any followup questions in your question and I will respond with an edit on my own.

I'm glad to hear you find my blog articles useful.

"How do you "heavy smooth" second harmonic when the base frequency could be 100x higher? "

Exponential smoothing has the property that the higher the frequency the more the tone will be attenuated.  If you are still getting multiple zero crossing then you still have higher harmonics present with significant amplitude.  Increase your smoothing parameter even more.  If you are working in floating point, you don't have to worry about destroying your fundamental frequency.  It will always be less dampened than the higher frequencies.

Since you have such a high sampling rate, you may also want to downsample your signal some as well.  There are many ways to do this, but since this is not an area of expertise for me, I will hold off on recommendations.  I can say the lowest calculation method is to just sample every nth point, but this does nothing to mitigate alias frequencies.  Another simple method is average intervals to get a new sequence.

The nice thing about picking a whole plus a half number of cycles for your frame is that the second harmonic (twice the fundamental frequency), and all even harmonics above it, are a whole (or very close to) number of cycles in the frame so they will be orthogonal to your basis vectors.  If there are a lot of harmonics present, you may want to increase the number of points above four per cycle to push the Nyquist frequency up for your frame length.

You can find my email address on my profile page.  Feel free to contact me there.

P.S.  Thanks to the upvoters.


How to Get Rid of Ripples from a Gradient Image of a Smoothed Image? I have a gray scale UINT8 image. When I smooth it with a Gaussian filter and then compute the gradient image with a Sobel filter in OpenCV, ripples would appear.

The first image is what happens when I use a 25x25 Gaussian filter with sigma value 4 and compute the second order derivative in the x direction with a size 3 Sobel filter.

The more I smooth the image, the more significant the ripples become.

The second image is what I got when I use a 49x49 Gaussian filter with sigma value 7. The Sobel filter is the same as the first image.

Why is the Gaussian filter creating more noise when it's suppose to get rid of them? How should I deal with this?

I think it happens due to 2 things:

Quantization
You are working using UINT8 Image, try convert it into floating Point Image.
You may do this by mO = im2double(mI) where mI is the UINT8 and mO is a floating point image in the range [0, 1]. You may also do it using mO = double(mI) / 255.
Gaussian Kernel Radius vs STD Ratio
Your Gaussian Filter Radius is 12 while $ \sigma = 4 $. This is a ratio of 3. Try higher ratios, something like kernelRadius = ceil(5 * kernelStd) this will make the Gaussian Filter result much smoother.


How to detect peaks ignoring soft declivity? I have to detect repetitions in a signal, and I am managing to do it by finding peaks and valleys, and cutting the desired pieces. The following figure shows the pattern wanted, shaded areas are the portions of interest:

The problem is: some of the signals that I'm analyzing, instead of starting with a constant zero, show a very smooth decline, as shown in the following figure:

As a consequence, a peak is detected as one of the first values, and the first repetition is detected wrongly. I'm looking for suggestions to avoid it. I have already tried (not sure if correctly):

Using derivatives (diff) to detect inclination.
Using a robust regression (robustfit) and finding inclination by the coefficients - too slow.

I'm using a MATLAB code created by someone else, modified by me (original can be found here):

function repeticoes = peakdet(v, delta, x)
%PEAKDET Detect peaks in a vector
%        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
%        maxima and minima ("peaks") in the vector V.
%        MAXTAB and MINTAB consists of two columns. Column 1
%        contains indices in V, and column 2 the found values.
%
%        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
%        in MAXTAB and MINTAB are replaced with the corresponding
%        X-values.
%
%        A point is considered a maximum peak if it has the maximal
%        value, and was preceded (to the left) by a value lower by
%        DELTA.
% Eli Billauer, 3.4.05 (Explicitly not copyrighted).
% This function is released to the public domain; Any use is allowed.

maxtab = [];
mintab = [];

v = v(:); % Just in case this wasn't a proper vector

if nargin < 3
  x = (1:length(v))';
else
  x = x(:);
  if length(v)~= length(x)
    error('Input vectors v and x must have same length');
  end
end

if (length(delta(:)))>1
  error('Input argument DELTA must be a scalar');
end

if delta <= 0
  error('Input argument DELTA must be positive');
end

mn = Inf; mx = -Inf;
mnpos = NaN; mxpos = NaN;

lookformax = 1;

for i=1:length(v)-1
  this = v(i);

  if isempty(maxtab)
    if this >= mx, mx = this; mxpos = x(i); end
  else
    if this >= mx, mx = this; mxpos = x(i); end
    if this < mn, mn = this; mnpos = x(i); end
  end

  if lookformax
    if this < mx-delta
      maxtab = [maxtab ; mxpos mx];
      mn = this; mnpos = x(i);
      lookformax = 0;
    end
  else
    if this > mn+delta
      mintab = [mintab ; mnpos mn];
      mx = this; mxpos = x(i);
      lookformax = 1;
    end
  end
end
    if size(maxtab, 1) == 0
        inicios = [];
        finais = [];

    elseif size(maxtab,1) >= size(mintab,1)
        inicios = maxtab(1:size(mintab,1),1);
        finais = mintab(:,1);
    elseif size(maxtab,1) < size(mintab,1)
         inicios = maxtab(:,1);
        finais = mintab(1:size(maxtab, 1),1);
    else
        disp('Erro ao reconhecer repeticoes!')
    end

    repeticoes = [inicios finais];

end

Any suggestions are welcome!

You mentioned that you have tried the signal derivative to catch the declines. What problems have you faced?

It seems that your problem might be addressed with different words: Let $s(t)$ be your signal, and let $u(t)=ds(t)/dt$. You want a detector to indicate when $a\le u(t)\le b$, with $a,b<0$, for at least $T$ seconds.

If this rewrite of the problem is true, the challenge now is to calculate $u(t)$. You can use u=diff(s), as you mentioned if the signal $s(t)$ has no noise, or you need to estimate $u(t)$ using other techniques like these two for instance.


DFT normalization for amplitude estimation After an interesting recent answer, I'm doing some research on proper DFT normalization for sinusoidal peak estimation. It's clear that to get the correct amplitude from DFT peaks we need to normalize by the window gain
\begin{equation}
S_w \triangleq \sum_{n=0}^{N-1} w[k]
\end{equation}
Which becomes $S_w = N$ for the rect window.

Now I'd like to find some rough demonstration, Harris has one (see equations 12 and 13 in the article below), but he chooses the easiest path by making exponentials cancel out. Well he actually does a pretty good job, it's me that thought it could be inferred to some more general relation.

For the rectangular window I've come up with this.
Let's say we have a sampled analytic sinusoidal signal $x[n]$, with frequency of $m/N$ cycles per sample (or $m$ cycles in $N$ samples), and compute the DFT $X[k]$

\begin{equation}
x[n] = A \ e^{j 2\pi m n/N}
\end{equation}

\begin{align}
X[k] & = A\sum_{n=0}^{N-1}e^{j \frac{2\pi n}{N} m} e^{-j \frac{2\pi
  n}{N} k} \\
\\
& =  A\sum_{n=0}^{N-1}e^{-j \frac{2\pi n}{N} (k-m)}\notag \\
\\
& = A \ \dfrac{1-e^{-j2\pi(k-m)}}{1-e^{-j\frac{2\pi}{N}(k-m)}} \\
\\
& = AN\delta \left[ k-m\right]
\end{align}

where $\delta[k]$ is the Kronecker delta.
To get the correct amplitude in freq space we need to divide by $N$.

Now, let's say we have a windowed signal, let me call it $x[n]$ again

\begin{equation}
x[n] = w[n] \ A \ e^{j 2\pi m n/N}
\end{equation}

and its transform, again overloading notation

\begin{equation}
X[k] = A\sum_{n=0}^{N-1} w[n] e^{-j \frac{2\pi n}{N} (k-m)}
\end{equation}

From here I'm not sure how to proceed, guess the result should be
$$
\require{cancel}\cancel{X[k] \overset{?}{=} A \ S_w \ \delta[k-m]}
$$

Edit: after some discussion with Robert I believe the previous equation is just wrong. It works for the rectangular window only, while for the other windows the only thing we can say is about the peak $k=m$, which can be easily demonstrated just from the sum.

  G. Heinzel et al - Spectrum and spectral density estimation by the
  DFT, including a comprehensive list of window functions and some new
  flat-top windows

  F. J. Harris - On the Use of Windows for Harmonic Analysis with the
  Discrete Fourier
  Transform

okay, since this is about scaling, i should be anal about definitions.

discrete-time signal: $x[n]$  where $n$ are only integer values.

DTFT:  $$ X\left( e^{j\omega} \right) \triangleq \sum\limits_{n=-\infty}^{+\infty} x[n] \ e^{-j \omega n} $$

it is necessarily the case that $ X\left( e^{j(\omega + 2 \pi)} \right) = X\left( e^{j\omega} \right) $ for all real $\omega$.

DFT:  $$ \begin{align} X[k] & \triangleq \sum\limits_{n=0}^{N-1} x[n] \ e^{-j 2 \pi \frac{k}{N} n} \\
& = X\left( e^{j\omega} \right) \Bigg|_{\omega=2 \pi \frac{k}{N}} \quad \quad \text{if } x[n] = 0 \text{ for } n<0 \text{ or } n \ge N\\ \end{align} $$

it is also necessarily the case that $X[k+N]=X[k]$ for all integer $k$.

and another fact is that $ X\left( e^{j0} \right) = X(1) = \sum\limits_{n=-\infty}^{+\infty} x[n] $

and $ X[0] = \sum\limits_{n=0}^{N-1} x[n] $

if x[n] is a "low-pass" function, we will expect $X[k]$ to peak at $k=0$.  and the value of that peak is shown as $X[0]$ above.

so, now consider a window function $w[n]$ that has DFT defined the same way

$$ W[k] = \sum\limits_{n=0}^{N-1} w[n] \ e^{-j 2 \pi \frac{k}{N} n} $$

normally we think of the window function as always non-negative, $w[n] \ge 0$, but it wouldn't necessarily have to be that.  nonetheless, just like $x[n]$, the window function has a DC component:

$$ W[0] = \sum\limits_{n=0}^{N-1} w[n] \triangleq S_w $$

and the DFT of the window function must also be periodic, $W[k+N] = W[k]$, for all integer $k$.

now there is a "translation property" of all forms of the Fourier Transform that has a particular expression for the DFT:

$$ \begin{align} \text{DFT} \left\{ x[n] e^{j 2 \pi \frac{m}{N}} \right\} & = \sum\limits_{n=0}^{N-1} \left( x[n] e^{j 2 \pi \frac{m}{N}} \right) \ e^{-j 2 \pi \frac{k}{N} n} \\
& = X[k-m] \\
\end{align} $$

for $m$ equal to an integer.

same applies for $w[n]$ with an amplitude scaling factor of $A$.

$$ \begin{align} \text{DFT} \left\{ w[n] \ A e^{j 2 \pi \frac{m}{N}} \right\} & = \sum\limits_{n=0}^{N-1} \left( w[n] \ A e^{j 2 \pi \frac{m}{N}} \right) \ e^{-j 2 \pi \frac{k}{N} n} \\
& = A \ W[k-m] \\
\end{align} $$

now we do expect the window function to be a "low-pass" function.  that's how we design them, so we also expect $W[k]$ to peak at $k=0$ (or any other integer multiple of $N$) and we expect the expression above to peak when $k-m=0$ (or $k=m$) and the value of that peak is

$$ A \ W[0] = A \ S_w $$

so if $A$ is the amplitude we're trying to find, we would take the value of the peak, wherever $k=m$, and divide that value by $S_w$ to get the amplitude of the complex sinusoid.

it's as simple as that, but it is also shown to be too simple.  we know that a real sinusoid is the sum of two complex sinusoid, so there will be overlapping interference (a.k.a. "leakage") of one onto the other.  we also expect that the frequency will not be exactly an integer number of cycles $m$ during the time of $N$ samples.  then the math gets a little more involved, but the scaling issues should, if $N$ is large enough, approach the result we have above.


Mathematical reason for MCPP in audio recognition? MCPP uses a weird graph of human phychoacoustics measurements called MEL scale to process digital audio.

Why is MEL used in DSP voice recognition if it's based on logarithmic human perception, and why is it effective?

MEL theory gives the ratio that average humans logarithmically hear equal distance pitches, and the cause for that perception ratio is unknown.

Is it because humans speak frequencies using the same frequency ratios by which they hear them, so that PC's need to know those ratios?

Is there a physical reason for the MEL to work well in DSP? i.e. is it an unwitting compromise between equal loudness contour relative to physical sound energy and the loudest most useful frequencies in voice recognition?

The Mel is an empirical model that describes the change in human pitch perception at different frequencies: until 1 kHz we perceive pitch changes in a linear way, after that, we perceive it in a logarithmic way, just as you said.

The pitch function relating real and perceived pitch differences is approximated based on the measurements made by Stevens in this work; an equation commonly used to approximate it is:

$$f_{mel}=2595\:\log_{10}\:10\:\left(1+\frac{f}{700}\right)$$

As I understand there is no mathematical reason for using it, is just choosing a model of how human beings perceive sound and applying it to some interactive human-machine application in order to enhance it's performance


Wavelet transform: How to compute the initial coefficients when only samples are available? In standard MRA we have that the space of functions at scale J can be expressed as
$$V_j = V_0\oplus \left(\bigoplus_{j=0}^{J} W_j\right)$$
where $V_0$ is spanned by the orthonormal system of the translates of the scaling function $\phi(x)$, i.e. $V_0 = \overline{\text{span}}\{\phi(x-k)\}_{k\in\mathbb{Z}}$ and the detail spaces $W_j$ are spanned by the wavelets ,i.e. $W_j = \overline{\text{span}}\{\psi_{j,k}\}_{k\in\mathbb{Z}}$, being $\psi_{j,k} = 2^{j/2}\psi(2^j x - k)$ and $\psi(x)$ the mother Wavelet. From the two scales equation and from the wavelet equation we obtain the iterative algorithm to compute the coefficients at the coarser scale from those of the finer one. That is, given a function $f(x)$ in $V_J$ (which can be for instance an approximation of another function up to such scale) then since by definition of MRA we have that $\{\phi_{j,k}\}_{k\in\mathbb{Z}}=\{2^{j/2}\phi(2^jx-k)\}_{k\in\mathbb{Z}}$ spans $V_j$, we can compute the coefficients $a_{J,k} = \langle f,\phi_{J,k} \rangle$ and then I can compute the coefficients at each scale $0\le j<J$ as
$$a_{j,k} = \sum_{l\in\mathbb{Z}} h(l-2k)a_{j+1,l},\qquad b_{j,k} = \sum_{l\in\mathbb Z}g(l-2k)a_{j+1,l}$$
where $b_{j,k} = \langle f,\psi_{j,k}\rangle$ and $h(n)$ and $g(n)$ are the scaling and wavelet filters, such that $$\phi(x) = \sqrt{2}\sum_{k\in\mathbb Z} h(k)\phi(2x-k),\qquad\psi(x) = \sqrt{2}\sum_{k\in\mathbb Z} g(k)\phi(2x-k)$$
So, once I have $a_{J,k}$ I can compute the whole decomposition, making filtering and all the operations I want and then re-synthetize the function inverting the iterative relation among the coefficients and rewriting
$$ f = \sum_{k\in\mathbb Z} a_{J,k}^{new} \psi_{J,k}$$
The problem is, in order to start the whole I need to compute the scalar products
$$a_{J,k} = \langle f,\phi_{J,k} \rangle$$
How I can I do it is I only have some samples of $f(x)$, say $\{f(x_i)\}_{0\le i\le N} $ ? (note that the samples can also be not equally spaced).
The solutions I heard about are
1) Use the low of large numbers and approximate
$$\hat{a}_{J,k} = \dfrac{1}{N}\sum_{i=1}^N f(x_i)\phi_{J,k}(x_i)$$
2) Solve the system of N equations
$$ f(x_i) = \sum_{k} \hat{a}_{J,k}\phi_{J,k}(x_i)$$
The first one is the simplest, and has the advantage that it can be computed very fast as a new sample arrives, however it has the problem that many coefficients will be undefined and many will be defined with very few terms. The second asks to solve a possibily big system of equations any time a new sample arrives. Are there any smart way to compute the tarting coefficients $a_{J,k}$?

I remember that using the discrete samples was termed  "wavelet crime" by Gilbert Strang (Strang and Nguyen, Wavelets and Filter banks, 1996, pages 232 sq.), or "wavelet sin" by others.

On prefiltering discrete signals to obtain "better" wavelet coefficients, here are some references I have in store:

Wavelet coefficient computation with optimal prefiltering, IEEE Trans. Signal Processing, 1994, Xia, X.-G. and Kuo, C.-C.J. and Zhang, Z.
On the initialization of the discrete wavelet transform algorithm, IEEE Signal Processing Letters, 1994, Abry and Flandrin
From the wavelet series to the discrete wavelet transform - the initialization, IEEE Trans. Signal Processing, 1996, Xiao-Ping Zhang and Li-Sheng Tian and Ying-Ning Peng
Sard-optimal prefilters for the fast wavelet transform, Numer. Algorithms, 1997, S. Ehrich
A new prefilter design for discrete multiwavelet transforms, IEEE Trans. Signal Processing, 1998, Xia, X.-G.
Quadrature prefilters for the discrete wavelet transform, IEEE Trans. Signal Processing, 2000, Johnson, B. R. and Kinsey, J. L.
Meaningful MRA initialization for discrete time series, Signal Processing, 2000, D. Veitch and M. S. Taqqu and P. Abry
Efficient Wavelet Prefilters with Optimal Time-Shifts, IEEE Trans. Signal Processing, 2003, Ericsson, S. and Grip, N.

It seems however that no prefiltering does not induce too much errors for compression, denoising, deconvolution. But for detection, or using multiple  different  wavelets at once (like with dual-tree wavelets, or unions of wavelet bases), or with multiwavelets, it is quite mandatory.


relationship between block size and detected frequencies in goertzel algorithm I'm trying to use the goertzel algorithm to detect a pair of frequencies in a signal.

The sample rate I'm using is 8000 hz, and the size of the incoming sample buffer is 60 samples.

I have a large array of 600 samples in size, so each time a new buffer comes in I move the data 60 samples to the right in the large buffer to make room for the new buffer, which I put at the front of the array (from index 0).

[ buffer0 | ... | buffer9 ]

I apply the goertzel algo to this large array. I'm looking for the frequencies 133.33 and 333.33 hz.

My algorithm gives me normalized magnitude. For the first frequency, I get a more or less stable detected magnitude around 0.5, which is also the amplitude of the original signal (I'm responsible for generating the signal myself). So that seems correct.

For the second frequency, I get a wildly varying normalized magnitude which seems to vary between 0.0 and 0.50, almost in a sine-like fashion. I know that the amplitude of this sine wave is 0.5.

When I save the generated signals and analyze them using audacity's spectrum feature I can see the correct frequencies in the spectrum, so it looks like the signal is being generated correctly.

My questions:

1) the algorithm seems to have a hard time detecting the second frequency. Why?
2) what is the relationship between the number of samples used for running the goertzel algo, and the detected frequencies?

My goal is to detect 64 unique combinations, so I was thinking of generating one frequency from a set of 16 and another one from a set of 4 and then just adding the two signals with 0.5 amplitude.

Thanks
B

Since Görtzel computes a part of the DFT, the relation between detectable frequencies is the same. So with 8000Hz sampling rate and 600 samples you cover the frequency range of -4000Hz to 4000Hz with 600 frequency intervals, and both your frequencies are points of the subdivision, thus both should be perfectly detectable.

Can you show the code for the Görtzel iterations?

Or even before that, did you build the test signals using $\sin(\pi\cdot f\cdot t)$ or $\sin(2\pi\cdot f\cdot t)$? Only the second gives the correct frequency, the first could explain your observations, since the 600 samples would then contain $12.5$ periods of the oscillation for $f=333.33Hz$, but full $5$ for $133.33Hz$.


Why do we use normalized angular frequency? Let's say we want to generate a signal in MATLAB, using sin, representing the A-major chord (440 Hz, 550 Hz, 660 Hz, each with amplitude 0.1).

I know we have to use the normalized angular frequency. But I don't understand why? Could someone kindly expalin this to me, please.

We use normalized frequency units in DSP mainly for convenience-- most DSP textbooks are written this way so that the theory can be presented in terms of essentially unitless quantities.

For instance, if you are talking about a sampled signal $x[n] = x(nT_s)$, where $T_s$ is the sample period, the discrete-time signal $x[n]$ now has no notion of time associated with it, because $n$ is an integer. Consequently, there can be no notion of "physical" frequency, in Hz or rad/s.

Practically speaking, what this means is that you end up having to "translate" frequencies back and forth, i.e. find out what cutoff frequency you want your filter to be in Hz, make sure your sample rate is greater than twice that, then find out what that frequency is in normalized units.

For your original question, you just need to make sure that your sample rate is high enough to capture all the frequency content. Your sample rate would need to be greater than $2\times 660 \textrm{ Hz}$ at a minimum, but practically you may want to make it even higher for a number of other reasons.

In MATLAB your script would look something like this:

clear
clc
close all

Fs = 2*700; % set sample rate higher if needed
Ts = 1/Fs;
t_duration = 1; % for now, one second of audio. can change if needed
t = 0:Ts:t_duration;
A = 0.1;
signal = A*sin(2*pi*440*t) + A*sin(2*pi*550*t) + A*sin(2*pi*660*t);

Note that you might want to change the phases of the individual sine terms to get the effect you are looking for.

EDIT 1:

Note that in the example code $t = nT_s = n/F_s$, therefore the argument to the sin function is:
$$
2\pi f t = \left( 2\pi\frac{ f}{F_s} \right) n
$$

EDIT 2:

Technically, in your example if all you are doing is generating a signal, you don't need to explicitly use normalized frequency; you just need to make sure that your sample rate is high enough. This answer was touching on the broader question of, in general, why someone would want to.


Applying DFT twice does not actually reverse an array. Instead, the first element stays in place while the rest of the array is reversed. Why? I've heard a million times that applying DFT twice will result in a reversed array, but that is not what actually happens. Instead, the first element remains where it was, and the rest of the elements are reversed. Is there an intuitive reason for that?
Octave code:
N = 16;
x = linspace(1, N, N);

fft(fft(x)) / N

Result:
1   16   15   14   13   12   11   10    9    8    7    6    5    4    3    2

The sequence is exactly what you should expect:
$$x[-n]=x[N-n]\tag{1}$$
Clearly, for $n=0$ $x[n]$ and $x[-n]$ have the same value.
It seems like you were expecting to see the sequence $x[N-1-n]$ instead of $x[N-n]$. If you wanted that sequence you would need to modulate the result of the first DFT:

N = 16;
x = 1:N;
c = exp(-1i*2*pi/N*(0:N-1));
x2 = real( fft( ( fft(x) .* c ) ) ) / N;    % real() just to remove rounding errors
x2 =

  16  15  14  13  12  11  10  9  8  7  6  5  4  3  2  1


Implementing QPSK Demodulation I conducted a research on how to implement a QPSK demodulator. I found out the Quadrature QPSK demodulation technique which seems very straight forward. Here is what i understood from the explanation of the block diagram, and would like to confirm from this platform if my understanding is right.  From the block diagram below, I am meant to do the following for a 2048 bits QPSK modulated signal:

Split the incoming 2048 bits of QPSK modulated signal into the I and Q phase
Multiply the 1024 bits in the I phase by the recovered carrier
multiply the other 1024 bits in the Q phase by the 90 shifted recovered carrier

My question is what do I do to this sets of I and Q phase outputs in order to get my demodulated bits ?

Most explanations I have seen just says it passes through a decision maker, but none has gone the extra length to explain what really happens within the decision maker

Your block diagram is correct. However, it seems to me that your understanding of how to implement it is not. There are no 'bits' in the modulated signal! The input signal is a noisy analog bandpass signal, and the I and Q components are derived from this signal as shown in the block diagram: by demodulating and filtering the incoming signal. Then this (complex) signal is sampled at the symbol rate (you may need timing recovery to know when to sample) and fed to a slicer. The slicer decides which of the four possible symbols was (probably) transmitted in each symbol interval. After the detected symbols are known you simply decode each symbol according to the coding in the transmitter (each symbol carries two bits).

EDIT: I didn't mention matched filters even though they were emphasized in Dilip's comment. For an optimum receiver they are of course necessary, but given the level of the question it seems to me that the block diagram is OK, and that the idea is probably that it should be implemented as is. Even though it will result in a sub-optimum receiver it can perform well under mild channel conditions. The only real difference is the type of low-pass filter used (matched to the transmitted pulses or just eliminating out-of-band noise).


Effect of Adding the cyclic prefix on the toeplitz matrix in OFDM Assuming we have $N$ symbols to transmit encoded in block $k$,

Performing $N$−iFFT at the transmitter, we now have

The resulted signal $x(k)$ has length of $N$.  inserting  a cyclic prefix $CP$ of size $D$, the length of signal will be $N+D$ instead of $N$.

Assuming we have channel $h$ of length $L$, the convolution of signal with channel can be written as:

$y = h*x_{CP}(k)$ = $Hx_{CP}(k)$ ,

where * is the convolution operation and $H$ is toeplitz matrix of size $(N+D+L),(N+D)$ built in matlab as below :

H = toeplitz([h(1) zeros(1,length(x_cp)-1) ], [h.' zeros(1,length(x_cp)-1) ]).';

As known, the signal $y$ has now the length of of $D+N+L$. However, the useful signal has the length of $N$ which is equivalent to $s(k)$

What I am asking about is the toeplitz matrix $H$ equivalent into $y$ after removing the delays $L$ and cyclic prefix $D$?  In other words, If I can write the $y$ in matlab as y = y(D+1:end-L+1); whose length becomes $N$ now, how can I write $H$ equivalent into this part ?

What you mean might be circulant matrix instead of toeplitz matrix. See section 3.4.4 in https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_chapter3.pdf about how the circular convolution in OFDM is represented by matrix operations (eq 3.130 onwards).

First, in almost all standard OFDM systems, you can assume $D \le L$. The cyclic prefix will be less than or equal to maximum multi-path delay, so as long as $D \le L$, the linear convolution $x * h$ gets converted circular convolution of $x$ and $h$.

When $y_c = H x_{cp}(k)$, you only need to take $N$ original elements for $x$.  $H$ is size $N \times N$. This is because you have already re-written the circular convolution in matrix form. Each row of $H$ will do dot-product with $x$ to generate $y_c[n]$. Like that there will be $N$ values of $y_c$ corresponding to each row of $H$.

$y$ is of length $N+L +D-1$. As you correctly mentioned $y_c = y(D+1:D+N)$.

Equivalent $H$ (size $N \times N$):

h(0)  0    0 ... h(L-1)  h(L-2) .. h(1)

h(1)  h(0) 0     0  ... h(L-1) h(L-2).. h(2)

h(2)  h(1)  h(0)  0 ..  0 0 h(L-1) h(L-2) .. h(3)

.

.

.

0 0 0 ...    h(L-1) h(L-2) .. h(2) h(1) h(0)

MATLAB command to generate above $H$:

Assuming $h$ is the vector representing channel having $L$ non-zero values.

H = toeplitz([h zeros(1,N-L)][h(1) zeros(1,N-L) flip(h(2:L))])

For a small example where $N=4, L=2, D=1$

$x = [x0, \,x1,\, x2,\, x3]$

$h = [h0\,,h1]$

$y = [x3h0\;, x3h1+x0h0\;, x0h1+x1h0\;, x1h1+x2h0\;, x2h1+x3h0\;, x3h1]$

$H$

h0 0 0 h1

h1 h0 0 0

0 h1 h0 0

0 0 h1 h0


Is there a "first" Fourier transform? Every time I've read "first Fourier transform", I've assumed it to be a misspelling of the fast Fourier transform. Recently though, I saw the term in both the title and keywords of an article from a seemingly serious, albeit not peer-reviewed, journal, Japanese Circulation Journal.

Development of a new analytical method for the electrocardiogram using short-time first Fourier transforms doesn't cite any resource or explain the algorithm. (fulltext PDF)

Is there a "first" Fourier transform?

Independently of what they call it, the process described in the paper is precisely the Fast Fourier Transform. This is made obvious by the FFT algorithm they describe, including the Gold-Rader bit-reversal procedure.

On a side note, if one tries Googling "SFFT", the S in the acronym stands for anything from sliding to spherical, but the F's always stand for Fast Fourier.


Implementing IIR filter using FSM How do I implement an Infinite Impulse Response (IIR) filter using Finite State Machine(FSM)? The problem is FSM has finite memory, but an IIR filter has infinite outputs. The closest match is this paper  . Unfortunately I do not have access to the paper.

Use feedback from output to input to perpetuate the output of the state machine.  You don't need an infinite number of states to produce an infinitely long output sequence.  For example, a state machine with a single memory element and feedback can be used to implement a first order IIR filter.  An example of such a system has the recursion relationship y[n] = 0.5y[n-1] + x[n] where y is the output x is the input and 0.5y[n-1] is the past output scaled, fedback and summed with the input x[n].


Convert decimal to Q5.27 Let's say I have a decimal like 0.9

To convert to Q0.27 I just multiply by 2^(27)

But with Q5.27, with 1 bit representing sign, how would I use the full range of bits available to me, to say convert 0.9 to a number between 16 and -16?

You can't really use Q5.27 in 32-bit storage. You need at least one bit for the sign, so I guess you'd call it Q4.27 since the sign bit is implicit.  The conversion from float is exactly the same as for Q0.27 except you can convert any numbers between -16 to just less than 16.

int(0.9 * 2^27)  = 0x07333333
int(-0.9 * 2^27) = 0xf8cccccd

As you can see you need the sign bit to represent a negative number.

-16 and 16 arise from:

0x7fffffff / (2^27) = 15.9999999925  (biggest +ve number)
0x80000000 / (2^27) = -16            (biggest -ve number)


How to properly time-shift a function in relation to another function I'm stuck on a time-variance question in my  homework and I think it has something to do with my understanding of how time shifting functions works.  The problem asks to determine whether the system defined below is time-variant or time-invariant.
$$ S: x(t) \to y(t) = x(t-2) + x(2-t) \tag{1}$$
I solved it such that
$$ x_{delayed}(t) = x(t - \tau)\tag{2}$$
and
$$ y_{delayed}(t) = x(t-\tau-2) +x(2-t +\tau)  \tag{3}$$

Then I delayed the original output $ y(t) $ from $(1)$ to
$$y(t-\tau) = x(t-\tau-2) + x(2-t+\tau) \tag{4}$$
which I believe is incorrect.  I think that instead, it should be
$$y(t-\tau) = x(t-\tau-2) + x(2-t-\tau) \tag{5}$$

I believe my problem comes from my understanding of how to delay a function related to another function.  Searching around didn't help much.  Ultimately my question is, should
$$ y(t) = x(-t) \tag{6}$$
$$ y(t - \tau) = x(-t + \tau\tag{7})$$
or
$$ y(t-\tau) = x(-t-\tau)\tag{8}$$
Using $(7)$ gives my my current answer, $(4)$, while using $(8)$ gives me what I believe is the correct answer, $(5)$.  I have found proofs for both answers which only confuses me more, so which is correct and how should I time shift the function? Any help is welcome and appreciated, thanks.

See that if you are going to go from $y(t)$ to $y(t-\tau)$ then you should replace every $t$ in the first one with $t-\tau$. This means that, in your case,

$$y(t-\tau) = x(t-\tau-2) + x(2-(t+\tau)) = x(t-\tau-2) + x(2-t+\tau)$$

In other words,

$$y(t) = x(-t)\implies y(t - \tau) = x(-(t - \tau))= x(-t + \tau)$$


Gaussian FIR filter with no multipliers? What is the simplest way to implement a Gaussian FIR filter with unity gain coefficients and no multipliers?

Please preface your answer with spoiler notation by typing the following two characters first ">!"

Note: A Gaussian FIR filter is an FIR filter with an impulse response that is a Gaussian function.  By "unity-gain coefficients" I mean all coefficients in the filter structure are 1.

This is an approximation, but you can make it as good as you like.

   Just use a cascade of several filters with rectangular impulse responses. In the simplest case this would be a two-tap filter. This works because of the central limit theorem. However, you will need to scale, because otherwise your resulting impulse response may become too large. The scaling could be done by bit-shifting.


what is the type of filter that has values ones? I used the filter with values ones but i don't know what kind of type is?

p=ones(1,L) where L is the length of the filter

in this case what we can name the filter?

This, on its own, would perform a summation operation of L elements of your input. If you modified it to p = ones(1,L)./L, then it would become a moving average (or in other words a low-pass) filter.

For more information please see this and this link.

Hope this helps


ROC of the function in the problem 9.14 of Oppenheim's Signals and Systems textbook I have solved the problem 9.14 in Oppenheim's Signals and Systems textbook, but my solution and the one in Slader is different. Problem is given above. And Slader solution is here.

I have also attached my solution below. Laplace transforms are the same but ROC in the Slader solution and mine is different. My question is how can this Laplace transform have an ROC as Re{s}>0. This should be wrong since this region contains the two poles of the Laplace transform which should not be the case for a rational Laplace transform.

Thanks in advance.

You're right, the ROC of the Laplace transform of a two-sided signal is a strip in the complex plane. In your case, the imaginary axis is inside the ROC, and the ROC is limited by the poles in the right and left half-planes. If the ROC were the right half-plane, the signal would be right-sided, which is clearly not the case.


Speech recognition using mfcc and neural network I am doing the project on speech recognition.I have implemented both the mfcc features and neural network.I am a bit confused in joining the two blocks.For example i have three different words with their '.wav' files.I extract 13 mfcc features from these three words respectively.There will be a matrix of mfcc e.g (13x100) for each word.I convert this matrix to a vetor.Than i will have three vectors of 1300 length (assuming each word is of same length).Than i pass these vectors as input to my neural network.I am confused about my training output.Whether it will be the same wav files from which i extracted the mfcc features or will it be something else.

Usually your training data and your test data are different, otherwise your results will be unrealistically good. However, for a first test of your implementation, you can simply use your training data as test material. In this case, the recognition performance should be very good. If this is not the case you know that something went wrong.


Implementing MVDR beamformer in the STFT domain I am trying to implement an MVDR beamformer for the first time. I was reading a couple of papers and books with many different notations and I am a bit confused.

In my case, without loss of generality, I have a spherical microphone array. Given that I produce the steering vector $d(f)$ with respect to a given 3D direction $\Omega_i=(\theta_i,\varphi_i)$, (azimuth and elevation). Let $y(f)=[Y_1(f), Y_2(f),\ldots,Y_M(f)]^T$ be the frequency representation of a single STFT frame. If I understood correctly my MVDR filter should be:

$$h_\text{MVDR}(f)=\frac{\Phi_y^{-1}d(f)}{\bar{d}(f)\Phi_y^{-1}d(f)}$$
Where, $\bar{(\cdot)}$ is the complex conjugate.

Now, my enhanced signal may be computed as:
$$Z(f)=h_\text{MVDR}^H(f)y(f)$$

I do this for each time frame, concatenate the result and follow to inverse STFT.

Is this a correct implementation?
If not, where was I wrong?
How do I estimate $\Phi_y^{-1}$ from my signal? Is that simply autocorr(y) in MATLAB?
Is there a nice python package with references or even one with an implementation for such a filter?

In this context $\Phi_y$ often describes the (estimated) power spectral density matrix  of $y(f)$, which is
$$\Phi_y(f) = E\{y(f) y^H(f)\} =
\begin{bmatrix}
 y_1 \cdot y_1^*(f),& y_1 \cdot y_2^*(f),& \dots    \\
 \vdots        & \ddots    &          \\
 y_M \cdot y_1^*(f),&               & y_M \cdot y_M^*(f)  \\
   \end{bmatrix},$$
where $M$ is the number of channels, and each element of the matrix corresponds to the cross-power spectral density of two channels for this frequency. Consistent estimates of PSDs in blockwise STFT processing can be obtained, for example, using Welch's method (pwelch in MATLAB).

I am also unsure about your steering vectors and complex conjugate transposes (which I assume is the $\cdot^H$ in your notation). They need to make the steering vector compensate the delay for a specific direction.

From my memory, I think the formulas should be

$$ h_\text{MVDR}(f) = \frac{\Phi_y^{-1}(f)d(f)}{d^H(f)\Phi_y^{-1}d(f)}$$

and

$$ Z(f) = h_\text{MVDR}(f)^T y(f).$$

I think the orginal MVDR paper is this one by Capon. You might want to double-check notations since a complex conjugate does change the beamformer's output quite a bit.

Edit: Thanks to Stanley Pawlukiewicz for pointing out my error!


What do the 1D filters represent when using imfilter? I am reading the source code of an algorithm that is used to process an image. While reading this source code (and others), I've found lines of code of the form

imfilter(image, [0.25 0 -0.25]', 'circular');

or

imfilter(image, [1 0 -1], 'circular');

I don't get what these kernels [0.25 0 -0.25]' (the transpose) or [1 0 -1] represent. Shouldn't kernels be 2D if the input is an image? What exactly do these specific kernels do to the image? I have seen several examples of imfilter being applied to an image with kernels of this form and I don't understand what the results of these operations should be (when I read the source code). Can someone provide some intuition?

If I'm not mistaken, a column vector will filter the image across its columns, treating each row independently of the others. Likewise, a row vector will filter across rows, treating all columns the same.

edit: Regarding an example - consider the simple image [1,1,1;0,0,0;-1,-1,-1]. It's constant along its rows (i.e., all the columns are the same) and a gradient along its columns. Let's filter it with a differentiating filter kernel [1,-1] along rows or columns. If we do imfilter(I,[1,-1],'circ') we obtain the zero image: since the image was constant along rows, filtering each row with a differentiating kernel gives the zero image. Oh the other hand, for imfilter(I,[1,-1]','circ') we obtain [1,1,1;1,1,1;-2,-2,-2]: each column gets differentiated independently and since all columns are the same, the resulting image is constant along rows.

As for the 0.25: this is merely a scaling of the whole image. You might as well filter with 1 as a filter weight and divide the result by 4, the effect is the same.


Relationship between the Real and Imaginary parts of a LTI causal system Prelude

I am writing an elaborate text on the relationship between the real and imaginary parts of a LTI causal system and how stability, causality and analyticity imposes various constraints on its magnitude response and the phase response. I am trying to explain three topics

Kramers Kronig Relationship
Bode's Gain-Phase Relationship
Hilbert Transform relationship

And finally, explain the connection between the three and their practical implications. This, I feel, is a much-needed text for people like me entering into this field. I have a list of questions, which I feel must come under the same context - please be patient in answering and kind enough to point out any mistakes.

Questions

The text, Kramers–Kronig Bode and the meaning of zero by John Bechhoefe, provides a short derivation of the gain-phase relation.

What is the motivation behind taking the logarithm of the response function in deriving the gain-phase relation?

Kramers Kronig (KK) relation requires only a causal system (analytic in the upper half of the complex $\omega$-plane) and does not talk about stability.

The Bode's Gain Phase relation requires a minimum phase system with no zeroes in the upper half of the complex $\omega$-plane and for non-minimum phase it becomes an inequality. In other words, Bode's Gain-Phase relationship necessitates that the response function must obey KK relations.
So, is the Gain-Phase relation a special case of KK relations accounting for stability too?
Causality in time domain implies analyticity in frequency domain, vice-versa.  Is my statement true?
How are Hilbert transform relationship and KK relationship related?
By Hilbert Transform relationship I mean the following (for a discrete case):
$$ H_{Im}(\omega)=-\frac{1}{2\pi}\int_{-\pi}^\pi H_{Re}(\lambda)\cot(\frac{\omega-\lambda}{2})d\lambda. $$
For a given response function $H(\omega)=H_{Re}(\omega) +j H_{Im}(\omega)$.
How do these relationships extend to higher dimensions? (for images and videos).

Thank you for your patience.

I believe the motivation is to make the math easier. Without the logarithm, you have $G(\omega)=|G(\omega)|e^{j\angle G(\omega)}$. By taking the logarithm, you have $\ln G(\omega)=\ln |G(\omega)|+j\angle G(\omega)$; by removing the exponential function and turning it into a sum, the relationship between $|G(\omega)|$ and $\angle G(\omega)$ becomes much easier to handle.
KK only requires causality and linearity. Bode is more strict because it also requires no zeros in the upper half of the plane (but make sure to look at note 16 in the paper).

As far as I can see, the value of the real part of the zeros is irrelevant; only their existance matters. So, stability is not really involved in this analysis at all.

The implication (causal -> analytic half plane) is true. I believe the "vice versa" (iff) is only true for stable systems, but I'm not 100% sure.

I don't know the answers to your questions 4 and 5.


Is there a technical term for this simple method of smoothing out a signal? Firstly, I am new to DSP and have no real education in it, but I am developing an audio visualization program and I am representing an FFT array as vertical bars as in a typical frequency spectrum visualization.

The problem I had was that the audio signal values changed too rapidly to produce a pleasing visual output if I just mapped the FFT values directly:

So I apply a simple function to the values in order to "smooth out" the result:

// pseudo-code
delta = fftValue - smoothedFftValue;
smoothedFftValue += delta * 0.2;
// 0.2 is arbitrary - the lower the number, the more "smoothing"

In other words, I am taking the current value and comparing it to the last, and then adding a fraction of that delta to the last value. The result looks like this:

So my question is:

Is this a well-established pattern or function for which a term already exsits? Is so, what is the term? I use "smoothing" above but I am aware that this means something very specific in DSP and may not be correct. Other than that it seemed maybe related to a volume envelope, but also not quite the same thing.
Are there better approaches or further study on solutions to this which I should look at?

Thanks for your time and apologies if this is a stupid question (reading other discussions here, I am aware that my knowledge is much lower than the average it seems).

What you've implemented is a single-pole lowpass filter, sometimes called a leaky integrator. Your signal has the difference equation:

$$
y[n] = 0.8 y[n-1] + 0.2 x[n]
$$

where $x[n]$ is the input (the unsmoothed bin value) and $y[n]$ is the smoothed bin value. This is a common way of implementing a simple, low-complexity lowpass filter. I've written about them several times before in previous answers; see [1] [2] [3].


Help identifying filter Apologies in advance for asking something so basic, this is not my field at all so I'm a bit lost where else to find this info. I'm trying to understand this piece of code that detects sudden drops in a signal, it goes something like this:
filter = [ones(1,10)*-0.05 ones(1,10)*0.05];

filterResponse = conv(exampleSignal, filter, "same");

Can anyone tell me what type of filter this is and why does it centers the signal at 0? Just anything to get me started onto more googling would be really helpful. Thanks in advance.

This is a moving difference filter which acts as a discrete time derivative and is less sensitive to high frequency noise. When all the samples in the waveform are the same (DC) the last ten samples will cancel out the first 10 samples for DC cancellation but when a change occurs within a 10 sample interval compared to 10 samples prior the difference will be maximized.
This filter is identical to a 10 sample moving average followed by a difference over 10 samples and provides a filtered estimate of the time derivative of the waveform as shown in the diagram below ($z^{-10}$ in the diagram refers to a 10 sample delay and MAF refers to a "Moving Average Filter" which is the sum of the prior 10 samples):

This equivalence form helps provide a more intuitive understanding of the filter with further details provided below:
The basic form of a differencing filter is $y[n] = x[n]-x[n-1]$, with coefficients given as [1 -1] and is a discrete time approximation of a derivative (as the Forward Euler mapping of $s$ which is the Laplace Transform of a time derivative):

Also consider how the time derivative is $\lim_{T \rightarrow 0}\frac{x(t+T)-x(T)}{T}$ in comparison to the form given above.
This filter has a simple high pass frequency response with the magnitude response shown below, given as the response from DC to the sampling rate where the frequency axis has been normalized by the sampling rate. Like the derivative which has a frequency response that increases as a function of f (given the Laplace transform of a time derivative is $s$), this filter is most sensitive to the highest frequency components in the signal (at the Nyquist frequency at half the sampling rate given by the normalized frequency of $0.5$):

If we insert 9 zeros to upsample the filter by 10, to get coefficients given as [1 0 0 0 0 0 0 0 0 0 -1], (which is the differencing stage in the first diagram provided above) such a zero insert creates periodicity in the frequency response, repeating the same frequency response ten times over the frequency range from DC to the sampling rate as shown in the magnitude response below (creating the classic "comb filter"):

The final response of the OP's filter is found from the convolution (or cascade) of the interpolated differencing filter with a 10 sample moving average filter (which removes the high frequency sensitivity as a low pass filter).
A moving average filter over 10 samples with coefficients given as [ones(1,10)] would have a Dirichlet response which is an aliased Sinc function. This is clear when we consider the coefficients of the filter is the impulse response of the filter, and the Fourier Transform of the impulse response is the frequency response. The coefficients are a sampled pulse, and the Fourier Transform of a pulse is a Sinc function. The response of a filter with coefficients [ones(1,10)] is shown below:

And the final form of the filter as the convolution of the coefficients of the interpolated differencing filter and the 10 sample moving average filter is shown with the response below (omitting the scaling by $0.05$). Convolution in time is multiplication in frequency, and we see that the resulting response is the product of the responses given above:


Calculating Savitsky-Golay filter coefficients - normalization Using the matrix equations found on Wikipedia I wrote some Python code to compute the Savitsky-Golay coefficients for $m$ data-points, using a polynomial of degree $k$. Here it is:

from numpy import *
from numpy.linalg import *

k = 4       # Degree of polynomial to use for regression
m = 25      # Window size, that is, the number of points. Must be an odd number.

z = range(-(m-1)/2, (m-1)/2+1)  # z is the distance from the central point.

jacobian, b = [], []
for i in z:
  for j in range(k+1):
    b.append(i**j)
  jacobian.append(b)
  b=[]

J = matrix(jacobian)

C = (J.T*J).I*J.T

I do realize that this code is in no way optimized and it could've been done in a much cleaner way, however I'm not a (professional) programmer nor a mathematician and would be really happy if I could get it to work at all.

My problem is that the coefficients computed become rather small, using the parameters above ($k=4$, $m=25$) on the order of $\pm 10^{-6}$. If I "go one step back" and multiply with the determinant I get really huge numbers ($\pm 10^{22}$).

I tried finding the greatest common divisor of these but got the result 2048 which leaves me with huge numbers (and it can't be right, result of numerical errors due to small/huge numbers or something?). I would like to have the coefficients in the form they're usually found in tables, with a normalization constant for each derivative, and an integer coefficient for each data point. Since I'm interested in the point of maximum curvature in a data series I want to calculate its 3rd derivative, I've been unable to find a table with these coefficients.

The SciPy library for scientific computing in Python contains functions for Savitzky-Golay filtering in its scipy.signal module since version 0.14.0, specifically you'll most likely be interested in the scipy.signal.savgol_coeffs function if you're only interested in the coefficients, or the scipy.signal.savgol_filter function which provides a nice interface for filtering alone.

When implementing some DSP method in Python, it's usually best to look into SciPy or scikits-* libraries rather than reinvent the wheel.


why does DFE feedback filter matters? I recently came across an article on the inventor of adaptative equalizer Robert Lucky, In this article : An Oral History by Robert Lucky (http://www.ieeeghn.org/wiki/index.php/Oral-History:Robert_Lucky).

And I have a couple of questions :

Quoting Lucky speaking about DFE :"It was the idea that you could use your own decisions to improve correction ability.  When you start transmitting data and the channel is very corrupted the bits you’re getting have a lot of mistakes in them. Nonetheless, you say, “Let’s assume these bits are correct and reconstruct what the perfect signal would have looked like with these bits. Compared to what the signal actually was, we’ll just try to make them look alike with our adaptive filter. ” You're basing this model on corrupted bits on the outside, so that the question is, “Is this ever going to converge ?” " ?

Can some one elaborate on this and expain how the DFE equalizer is assure to converge and provide better equalization than the simple forward filter system ?

My second question is about the design of the filter, I choose to use an MMSE criterion, is that the best alternative ?
My DFE equalizer is for HF ionopheric channels, how can I update the filters coefficients, i.e., based on the channel estimate, what is the exact algorithm ?

Any Hints, thoughts, or ideas will be very appreciated

I can't provide a lot of theory but I can provide this qualitative description of the problem:

Normally the output of your digital filter sampled every symbol period gives a point on your I/Q plane. Decision feedback usually selects the constellation point closest to this point and modifies your filter to make the output(given the same input) closer to the given constellation point.

If you select the correct point each time then your filter will hopefully converge(approximately) to the set of coefficients which minimize the above distance.

There are some obvious issues with this..

if you aren't selecting the correct constellation point then your filter is converging to the wrong thing. More likely(lets say in a hypothetical AWGN) you will select a point sometimes right, sometimes wrong - hopefully if you select it right more than wrong, the rights will dominate the filter updates.
Your transmitted pulses need to be mostly evenly distributed among the constellation, otherwise I think you could have a DFE converge to something perfect for one pulse train but maybe terrible for the others.

Intuitively you're attempting to invert the frequency response of the channel or reverse the impulse response but this is impossible for finite length equalizers.

I don't know what is currently used but I'm not sure if the standards which use OFDM have DFE or just do channel estimation and lock it down until the next pulse.

Your question about convergence: I don't know, I'm not sure how to analyze DFE, the adaptive filter course I took didn't address it - I have my doubts on how easy a proof is showing it will converge because even proving relatively simple things in adaptive filter theory is really hard(like the convergence pro0f of LMS algorithm was nontrivial even though it's a super simple filter).

I also have seen multiple proofs of algorithms like CMA which either proves it does converge or doesn't converge, or converges for a while and then diverges with a set of assumptions I'm not sure make the proofs relevant to real world applications.

Here's a short MS thesis which talks about some of this:

https://pdfs.semanticscholar.org/106d/7871c03d837a57be6cdf755d4efd44f5e96c.pdf

Basically, bad things can happen.


How to remove dark shadows from camera capturing when combining LED spots and ambient light I'm getting strange moving shadow strips on my camera capture during a computer vision experiment:

As the gif image above shows, those strips move to bottom slowly. I have also noticed that it happens only if I use a combination of ambient light and two 10 Watt LED spots aside the camera. Using only ambient light or only LED lights or ambient light plus a single LED spot the problem doesn't happen:

I'm using OpenCV to grab 320x240 images from the camera at 15 fps. The camera model is a Microsoft LifeCam Studio Q2F-00013. The 2 LED spots are like this one:

and power source is 220 volts at 60Hz.
I'm looking for a way to get rid of those shadows. Using v4l2-ctl I can see some camera controls to tune the camera:
$ v4l2-ctl -d 2 -l
                     brightness 0x00980900 (int)    : min=30 max=255 step=1 default=133 value=133
                       contrast 0x00980901 (int)    : min=0 max=10 step=1 default=5 value=5
                     saturation 0x00980902 (int)    : min=0 max=200 step=1 default=103 value=103
 white_balance_temperature_auto 0x0098090c (bool)   : default=1 value=1
           power_line_frequency 0x00980918 (menu)   : min=0 max=2 default=2 value=2
      white_balance_temperature 0x0098091a (int)    : min=2500 max=10000 step=1 default=4500 value=4500 flags=inactive
                      sharpness 0x0098091b (int)    : min=0 max=50 step=1 default=25 value=25
         backlight_compensation 0x0098091c (int)    : min=0 max=10 step=1 default=0 value=0
                  exposure_auto 0x009a0901 (menu)   : min=0 max=3 default=3 value=3
              exposure_absolute 0x009a0902 (int)    : min=1 max=10000 step=1 default=156 value=156 flags=inactive
               pan_absolute 0x009a0908 (int)    : min=-529200 max=529200 step=3600 default=0 value=0
              tilt_absolute 0x009a0909 (int)    : min=-432000 max=432000 step=3600 default=0 value=0
             focus_absolute 0x009a090a (int)    : min=0 max=40 step=1 default=0 value=8
                 focus_auto 0x009a090c (bool)   : default=1 value=0
              zoom_absolute 0x009a090d (int)    : min=0 max=317 step=1 default=0 value=0

But sincerily I don't know where to start to tweak it in order to remove those shadows.
EDIT:
According to this power_line_frequency=2 means 60 Hz which indeed is the power line frequency here. Playing with exposure_auto/exposure_absolute also don't provide any improvement. The only setting which worked so far have been setting brightness to 173. Anyway, I don't know if increasing the brightness is actually a robust solution.

Either the problem is simple and the powerline frequency is just set wrong. There is also a chance that the lights use PWM to drive the LEDs so they basically blink at high frequency that beats with the rolling shutter rate, so it is possible that those lights are unsuitable for cameras with rolling shutter.


Can I use frequency hopping spread spectrum technique for audio steganography? I have been trying to use spread spectrum technique for audio steganography. I read a lot about it.But i didn't find any information about the steganography scheme using FHSS(frequency hopping spread spectrum). Isn't it possible to do steganography in .wav files using FHSS?

Yes, it is possible to hide data inside audio files using FSSS technique. Have a look at this publication and this one.

Since HAS (Human Auditory System) is more sensitive than HVS (Human Visual System), you are usually limited with audio steganography comparing to image steganography. There are not many detalied publications with "Audio Steganography" tag. Interests about hiding digital data inside audio files have tended into watermarking field. So I strongly suggest you to research whatever you are looking for as "Audio Watermarking". You will find more contents related to watermarking.

Few more words about Audio Steganography:

LSB technique provides you more capacity to hide inside audio files, yet it will be very fragile to resist against any attacks such as mp3 compression. According to my experiences, DSSS and Echo Hiding methods are more robust comparing to LSB Coding and Parity Coding. DSSS provides a bit more capacity than Echo Hiding, but the noise is audible comparing to Echo Hiding.


What is the distribution of it? If $\theta$ is uniformly distributed in $(0, 2\pi),$ then what is the distribution of $e^{i\theta},$ where $i = \sqrt{-1}?$ And what are the statistical properties of $\left[e^{i0\theta}\, e^{i1\theta}\, e^{i2\theta}\dots\, e^{i(N-1)\theta}\right],$ where $\theta$ is a single observation and $N$ some integer?

$e^{ik\theta}$, where $k$ is any non-zero integer, is uniformly distributed on the complex plane unit circle. The expected value of $e^{ik\theta}$ is $\displaystyle\frac{\int_0^{2\pi} e^{ik\theta} d\theta}{2\pi} = \frac{i - ie^{2\pi i k}}{k} = \frac{i - i}{k} = 0.$ The variance of $e^{ik\theta}$ is $1,$ because all values that may be observed are at a distance of $1$ from the expected value $0.$

$e^{i0\theta} = 1$ is constant.

If you by notation $(0, 2\pi)$ mean that values $0$ and $2\pi$ are never observed, then one may observe $e^{ik\theta} = 1$ only if $|k| > 1$, because the distribution wraps around passing the point $e^{ik\theta} = 1$ at $\theta = 2\pi/k.$


Assumptions for Hurst exponent calculation Are there any general assumptions for the calculation of the Hurst exponent?

Does the signal need to be stationary, for example?
Does it depend on the method?
What about the length of the time series, is longer better?

I am interested in using implementing it for the analysis of an EEG signal. I have read many articles but can't find answers for these questions.

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

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


Describe quality of rectangular signal detection I have some process, which produces some rectangular signal. I have an algorithm, which recovers this signal. To test my algorithm I compare the recovered signal with known:

My question is:
Could You help me with some common and well-established methods to describe the quality of my detector?
Of course, I am able to calculate the number of non-detected jumps, delay its deviation. The question is that I don't want to use homemade solutions. Hope there should be a big science with well-established methods and criteria to perform such tests.

It is a rare individual who develops detection algorithms as a recreational pursuit.  If you are not such a person then the reason for your algorithm should provide, what is favorable and not favorable in some measurable way.

Your question has some ambiguity because random delay, missed transitions, and false transitions is all you have mentioned. "Big Science Tests" are based on probabilistic assumptions (a model) that applies.  If your transitions were random with some rate $\mu$, a Poisson model would be promising.  One might look at why you don't see a one when there is a one which is the error mechanism, which typically in Signal Processing is associated with noise.

Detection is usually cast as a Hypothesis Test with alternative hypothesis $\mathbf{H_0}$ and $\mathbf{H_1}$  where one hypothesis is signal present and the other is signal absent.

Another sort of hypothesis can be tied to a parameter. As an example, A Poisson Process has a mean rate $\mu$ which can be $\mathbf{H_0}$ and $\mathbf{H_1}$ can be the case the observed mean rate $\mu_{observe} \; \ne \;\mu$

As Mbaz suggests in his comment, missed and false detections can be combined as an error probability which might look like

$$P_{error}= P( \text{ones received}| \text{zeros sent}) P(\text{zeros sent})+ P( \text{zeros received}| \text{ones sent}) P(\text{ones sent})
$$
if the ones and zeros are independent.

The essentail point is that you need to understand your data model to use a standard test.

This is not an exhaustive list of types of statistical tests you can use. and many in signal processing have intractable forms. There are paired tests. There are non parametric tests. There are robust tests. Statisticians are inventing new tests every day.

Your data model on the other hand, is up to you.

Lastly, making things up is actually common and as long as it isn't fake, is often acceptable.


Phase and Amplitude estimation Context:

I have used FFT many times, but for real, non-periodic signals I consider it a poor estimator.

For most of my applications I am only interested on the power spectrum, so I use the Welch's method. I really like its approach to take overlapping segments, apply a window to each of them, calculate the FFT and finally average the result to get the power spectrum estimation. The method is easy to understand and its results are smooth and stable.

However I have one application that I am also interested in the phase of the signal. In fact am not interested in the phase estimation of particular signal segment, I want to see how the phases of very similar spread spectrum signals vary over time.

My signals spectrum are spread, but they always have a peak in the PSD, however this peak frequency varies over time. I will probably only need to trend the phase related to the frequency of this peak.

Questions:

1) Does it make sense to use Welch's method for magnitude and phase, instead of power? I mean average magnitude and phase, instead of power.

2) Is there any magnitude and phase spectrum estimator that is smooth/stable like Welch's method?

answer to 1, No,  the  phase of an tone isn't constant from one FFT to the next, so averaging phase provides no reduction of error unless there is reason to think there is like in a coherent algorithm like SAR.

Answer to 2 is a bit more complicated but is yes, because you need to do carrier and timing recovery in a communications systems

For something less structured There are several frequency estimation algorithms that work on FFT data like:

  Macleod, Malcolm D. "Fast nearly ML estimation of the parameters of
  real or complex single tones or resolved multiple tones." IEEE
  Transactions on Signal processing 46.1 (1998): 141-148.

and this just a personal favorite. There are 270 papers that cite that paper and many before that.  This method is nearly ML, so not Welch but not shabby.

  La Scala, Barbara F., and Robert R. Bitmead. "Design of an extended
  Kalman filter frequency tracker." IEEE Transactions on Signal
  Processing 44.3 (1996): 739-742.

Is more along the lines of tracking frequency and there are 170 papers that cite this paper.

a good book is:

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

A recent tutorial paper

Vila-Valls, Jordi, et al. "Are PLLs dead? A tutorial on Kalman filter-based techniques for digital carrier synchronization." IEEE Aerospace and Electronic Systems Magazine 32.7 (2017): 28-45.

might be useful.

Phase introduces complications because phase is relative to some reference and having the same reference for multiple tones may or may not be meaningful.

One of the regular contributors here, Peter K has a web page (which I don't recall where) that has a lot of material that may be of interest.


Difference between 2$\pi f$ and $\omega$ in Fourier transform What is the difference when we use $e^{-j2\pi f}$ and $e^{-j\omega n}$ for Fourier transformation?

You might be mixing two concepts, pertaining to (following Robert Bristow-Johnson)

the "analog context" (more formally, the "continuous-time case"),
the "digital context" (more formally the "discrete-time case").

The first concept corresponds to the continuous Fourier transform, for which you can use
a form of  normalized frequency cycles per second or Hertz ($e^{-j2\pi f}$), while the second is the angular frequency in radians per second ($e^{-j\omega}$).

The second one, with integer index $n$ relates first instance to discrete summations, e.g. the discrete-time Fourier transform (DTFT), where the Fourier kernel writes  $e^{-j2\pi f n}$ with $f\in \mathbb{R}$ (for instance) or $e^{-j\omega n}$ with $\omega\in \mathbb{R}$ (for instance).

The first concept is in use in the inverse DTFT, since it involves integrals.


concatenating short noise sequences I'm trying to merge 2 identical audio files. The Audio file is nothing more than white noise generated using MATLAB commands

block_len = 4096;
Y = wgn(block_len, 1, 5, 120, 'dBm', 'real');

Rather than generating a longer block of data, i would like to make a longer file by concatenating the block of 4096 samples.

Therefore I tried a simple,

datat = [Y;Y;Y;Y;Y];
soundsc(datat)

The output file is not smooth enough. There is an audible periodic pattern that indicates the multiple 'Y' segments. I believed this could be due to the imperfect combination at the beginning and End of each Y segment.

Therefore I added the lines

Y(1) = 0;
Y(end) = 0;

Concatenating after these still does not have a smoothness to it. I searched on stackoverflow and found this link.

I tried the code based on the link above. I still can't get a smoothness to the white noise. What other methods can i try to make the sound as smooth as possible?

Rather than generating a longer block of data, i would like to make a longer file by concatenating the block of 4096 samples.

Bad idea. That means your noise becomes perfectly correlated with a period of 4096 samples, and that's definitely not white noise anymore, and you'll stand a realistic chance of noticing that audibly; depending on the sampling rate, that $\frac{f_\text{sample}}{4096}$ might be a directly audible frequency.

A lot of signal processing algorithms are sensitive to periodicity, so: really, don't do that.

By the way, you might be clipping. Gaussian noise amplitude is unbound, but your sound card's acceptable numeric range very much is.

  The output file is not smooth enough. There is an audible periodic pattern that indicates the multiple 'Y' segments.

Which was predictable, because that's exactly what you've been asking it to generate.

  I believed this could be due to the imperfect combination at the beginning and End of each Y segment.

No, this was you; you introduced periodicity, and your human sound system notices that. Congratulations! Your hearing is functional :)

You using a "blend over" technique on white noise says you  haven't really thought about what white means: it means that there's no correlation between successive samples. Thus, if you concatenate two different, uncorrelated white noise sequences, you get a white noise sequence.

If you concatenate the same sequence, you get something periodic.

What you need to do is  generate white noise of the desired length instead of concatenating short sequences. You forgot to mention why you would even want to do that, because (haven't benchmarked matlab) white noise generation should really be fast.


How to choose the right segment length for a given signal? Im working on plotting the power spectral density of a random time series measured. I know the target spectrum.The idea is check the measurement vis-a-vis the target spectrum. However, I get slightly different plots for segment lengths of 128 & 256, its noisy for higher segment lengths. For 128 its a closer match. For 256 its smooth, but there is a slight bulge in the spectrum at the top (when compared to the target).

Im not sure which one to go with to validate the measuring device.
Please advise.

Edit : 1. Hand sketch of the look of the plot added for 256 point segment length.
       2. For 128 point, its more or less an exact match.

The DFT in of itself is not a consistent estimator, or in other words, the variance of your transform doesn't decrease as the transform length increases.

Another way to think of it is that the DFT is a linear invertible transform, so whatever "information" is contained in the time domain is contained in the frequency domain.  The segment that is 256 points is expected to be less smooth than a 128 point segment. If the longer segment has more variability, it's DFT will also have more variability.

To have a consistent estimator one can use the technique attributed to Welch, which means that you average a a set of DFT segments and that generally means averaging bin magnitudes (or magnitude squares).

Picking a correct segment length is related to how many segments you will average.  The number of averages is often tied to the amount of data you can average, as well as a target bin variance.

If the data has a spectrum that varies, like a speech signal, one can average too much and obscure the dynamic nature of the data.

Window functions have a large influence on smoothness, and are implicitly related to segment length.

Without knowing what your target spectrum's purpose is, any notion of  "correct length" will be subjective.   Subjectivity isn't necessarily bad.


Proper HW circuit to decode PWM signal using Arduino Mini I want to properly connect a Radio Control Receiver signal to a Arduino Pro Mini, Pin 11 so I can detect the 1 ms to 2 ms PWM signal.  I calculate a value of 220 ohms will limit the circuit current to 22mA.  The Arduino strongly suggests using 20mA with a 40mA max on input pins.  Do I need a .1uF cap to ground?  Electrolytic cap?
 I would think adding a cap will cause the signal line to ride at a certain voltage level instead of having a sharp leading/trailing signal levels.
TIA

Discretize Amplitude or Time – that is the (signal processing question)

Seeing that you're using a microcontroller running at at leas 8 MHz (probably 16 MHz): It's way easier to measure the time a signal was high than to first average out the PWM  cycle in analog domain and then convert that analog voltage to digital.

If this was a sigma-delta ADC, which are pretty popular for most applications these days: Technologically, the conversion from analog voltages to digital numbers happens through counting the time a signal was high; I think the Atmega328p uses a successive approximation ADC, though.

The technological reality is that for signals below a couple of GHz, it's easier to build fast digital silicon logic than to build exact analog silicon. Thus, the less complex, more accurate and often less power-hungry solution is always to rather measure in time direction than in amplitude direction.

Implementation

So, I'd stay with the 220 Ω resistor: It really limits the current in the misconfiguration case, and is negligible compared to the proper input impedance.

Whether or not to add a capacitor really depends on what amount of noise you're expecting: you'll need a lot of noise to convert a stable TTL "high" to a "low", or vice versa. It still makes sense to do a bit of filtering, but not really with 0.1 µF:

----|220 Ω| ----·------ GPIO PIN
       R        |
                = C?
                |
                v GND

If we understand the RC combination as low-pass filter, then 220 Ω -- 100 nF leads to a cutoff frequency of ca 8 kHz, so not enough to average out a PWM cycle (if you wanted to go the analog voltage route), but maybe a little too much if you wanted to the low/high edges nice and sharp for detection. Assuming no significant inteferers between that frequency and 80 kHz, I'd maybe start with 10 nF – I think the ATmega328p has smitt trigger inputs, so you're really pretty safe from transients.

So, you only need to implement a way of measuring time: for the large time scales you're talking about, this could easily done by waiting for the pin to become high, and then just counting the CPU cycles until it turns low again. But that blocks the CPU for anything else – it'd be more elegant to use to use the built-in Timer/Counter units of the microcontroller, which have been invented for exactly that reason. You can program them to start counting at the rising edge of an external input signal, and stop on a negative edge.


Downsampling a signal with decimate I'm experimenting with decimating a signal, in this case a unit impulse.

I'm using Python, with pylab.  First, I create a unit impulse, and decimate it by 5.

x = r_[zeros(0), 1, zeros(100)]
N = 2 ** 14
q = 5

y = decimate(x, q, ftype="fir")
subplot(211)
title("Original")
stem(range(len(x)), x)
subplot(212)
title("Decimated - FIR")
stem(range(len(y)), y)

figure()
subplot(211)
semilogx(log(abs(fft(x, N))))
subplot(212)
y = decimate(x, q, ftype="fir")
semilogx(log(abs(fft(y, N))))

This results with the following plots

I then add a few samples of delay before the impulse, by changing x to:

x = r_[zeros(3), 1, zeros(100)]

This results in the following plots

In the second set of plots, the resulting decimated signal is no longer a single sample, but has been distorted.

If I delay the signal with 5 - and any multiple of q - samples, I get the first set of plots again.

The source code for the decimate function is, https://github.com/scipy/scipy/blob/master/scipy/signal/signaltools.py#L1570

def decimate(x, q, n=None, ftype='iir', axis=-1):
    if not isinstance(q, int):
        raise TypeError("q must be an integer")

    if n is None:
        if ftype == 'fir':
            n = 30
        else:
            n = 8

    if ftype == 'fir':
        b = firwin(n + 1, 1. / q, window='hamming')
        a = 1.
    else:
        b, a = cheby1(n, 0.05, 0.8 / q)

    y = lfilter(b, a, x, axis=axis)

    sl = [slice(None)] * y.ndim
    sl[axis] = slice(None, None, q)
    return y[sl]

I'm using a fir low pass filter before decimating, the impulse response of the filter is

This explains why the impulse is distorted when there is a delay, the decimation is selecting parts of the impulse response, when the delay is a multiple of the decimation, it only selects the zero's of the impulse response, and one non-zero sample at the peak.

Is there a way to decimate a unit sample with an arbitrary delay, which results in a scaled unit sample output?

You seem to understand correctly what is going on. However, I'm not sure what you expect to get. Take your initial example. Let your input signal be $x[n]$:

$$
x[n] = \delta[n]
$$

The first step in the decimation process is that the input signal is convolved with the anti-aliasing filter's impulse response $h[n]$:

$$
\begin{align}
x_f[n] &= x[n] * h[n] \\
&= \delta[n] * h[n] \\
&= h[n]
\end{align}
$$

Next, the filtered signal is downsampled by a factor of $q$ ($5$ in your example).

$$
x_d[n] = x_f[qn] = h[qn]
$$

As you noted, for FIR filters whose orders are a multiple of $q$ (actually, since the filter is linear-phase, the order need only be a multiple of $\frac{q}{2}$), the taps at time delays $qn$ are zero for all $q \ne 0$. Therefore, $x_d[n]$ is nonzero only at $n=0$, as you found.

When you apply a time delay to the input impulse $x[n]$, then the filtered output is just delayed by the same amount, since the filter is linear and time invariant (LTI):

$$
x_f[n] = h[n-D]
$$

$$
x_d[n] = x_f[qn] = h[qn-D]
$$

Again, as you noted, this has the effect of plucking out a different set of taps from the filter's response, such that the decimated output signal is no longer zero for all but one sample (i.e. it doesn't look like an impulse any more). This is to be expected. Why?

Remember that the discrete-time Fourier transform (DTFT) of a discrete impulse is flat (in magnitude) in the frequency domain. If $x_d[n]$ is to be equivalent to a delayed impulse, it must also have flat magnitude in the frequency domain. However, its DTFT just a scaled copy of the filter's frequency response:

$$
x_d[n] = h[qn-D] \Leftrightarrow e^{-j\omega D} H\left(\frac{\omega}{q}\right)
$$

where $H(\omega)$ is the filter's frequency response. In order for the decimated output $x_d[n]$ to be equal to a delayed impulse, the filter must have a brick-wall response that is perfectly flat across its passband (up to the post-decimation Nyquist frequency) and zero everywhere else (so no aliasing leaks back in after downsampling, making the result non-flat). This isn't realizable unless you have an infinite amount of time and resources.

Since the filter is the source of the "distortion" that you don't want, you might consider trying the process again without a filter. But, consider what you would get then:

$$
x_f[n] = x[n] = \delta[n-D]
$$

$$
x_d[n] = x_f[qn] = \delta[qn-D]
$$

If $q$ is not a multiple of $D$, then $x_d[n] = 0\ \forall\ n$, which is probably not what you would want either.


MATLAB FFT Amplitude I am trying to understand FFT, DFT through Matlab and I am fairly inexperienced in all of these. I have a time vector t and a corresponding simple sine wave x = 0.25*sin(t). Now I would like to view the 1-sided FFT of this

%MATLAB
%define time and signal
delta_t = 2.6;
t = 0:delta_t:600;
x=0.25*sin(t);
L=length(x);

%define frequency vector
f_low  =1/t(end); %lowest frequency for this range
fs = 1/delta_t; %sampling frequency
f=[1:L]*f_low;

%calculate FFT
FFT = fft(x);

%double sided spectrum
P1=abs(FFT)/L;

%single sided spectrum
P2=P1(1:L/2+1);
P2(2:end-1)=2*P2(2:end-1);
f2=f(1:L/2+1);

%plot
figure()
plot(f2,P2,'LineWidth',1.25);
grid on
xlabel('f')
ylabel('abs(FFT)')

The result of this varies with delta_t (all the while respecting Nyquist theorem). The amplitude of the FFT is different for differnt delta_t and ideally shouldn't it be 0.25, equal to the amplitude of the wave?

Can someone please explain me why does the amplitude vary with delta_t?

Thank you!

Your code has some troubles, and perhaps also your theoretical understanding. Let me put here a mini summary of observing a practical sine wave and its frequency spectrum using FFT function of MATLAB / Octave.

Assume that there's a continuous-time ideal infinitely long sinusoidal wave with frequency $\Omega_0$ given as:

$$ x(t) = A \cos(\Omega_0 t) , \tag{1} $$

Let's sample this infinetely long signal within a finite duration of $[0 , T_s]$ using a sampling rate of $F_s$ or the sampling period of $T_s = \frac{1}{F_s}$, so that $t = nT_s$, and  we obtain $L$ samples of the signal :

$$ x_L[n] = x[n]w[n] = A \cos(\omega_0 n) ~~,~~n = 0,1,...,L-1 \tag{2}$$

where $\omega_0 = \Omega_0 T_s$ and $w[n]$ is a rectangular window, used to represent the finite duration of observation.

$$ w[n] = \begin{cases} { 1 ~~~,~~~~ 0 \leq n \leq L-1 \\ 0 ~~~,~~~ ow }\end{cases} \tag{3}$$

DTFT (Discrete-time Fourier transform) modulation property can be used to analyse the frequency domain representation of $x_L[n]$ :

$$
X_L(e^{j\omega}) = \frac{1}{2\pi} X(e^{j\omega}) \star W(e^{j\omega})  \tag{4}
$$

where $W(e^{j\omega})$ is the DTFT of the rectangular window :

$$
W(e^{j\omega}) = e^{-j \omega (\frac{L-1}{2}) } \cdot  \frac{ \sin(\omega L /2)}{\sin(\omega/2)} \tag{5}
$$

A plot of this function is shown below :

Again we also we know that the DTFT of $x[n] = A\cos(\omega_0 n)$ is :
$$
X(e^{j\omega}) = A\pi \delta(\omega + \omega_0) + A\pi \delta(\omega - \omega_0)  \tag{6}
$$

Hence we see that the DTFT of the finite length observation $x_L[n]$ is

$$
X_L(e^{j\omega}) = \frac{A}{2} [ W(e^{j(\omega+\omega_0)}) +  W(e^{j(\omega-\omega_0)})] \tag{7}
$$

Now, when you take $M$-point FFT of the $L$-point sequence $x_L[n]$, then you are effectively calculating the $M$-point DFT (discrete Fourier transform) which is equal to $M$ uniform samples of DTFT $X_L(e^{j\omega})$  as given by

$$
X_L^M[k] = X_L(e^{j\omega})|_{w = \frac{2\pi}{M} k} ~~~, ~~~k = 0,1,...,M-1 \tag{8}
$$

which yields the complex frequency vector that you compute using FFT function:

$$
X_L^M[k] = \frac{A}{2} [ W(e^{j(\omega_k+\omega_0)}) +  W(e^{j(\omega_k-\omega_0)})]  ~~,~~k = 0,...,M-1 \tag{9}
$$

So with an M-point FFT of $x_L[n]$ we calculate the samples of the shifted versions of $W(e^{j\omega})$.

An example plot of $|X_L^M[k]|$ for $L =  16$ , $M = 1024$, $A = 5$ and $\omega_0 = 0.2 \pi$ is shown.

Note that the peak values (which are $A \cdot L /2$ when $\omega_0$ meets some criterion, see below) are obtained for some $k=k_0$ and $k = M-k_0$. Consider the first peak (to the right of the plot here) and it can be shown that its value is

$$ |X_L^M[k_0]| = |0.5 A [  W(  \frac{2\pi}{M} k_0 - w_0) +  W(  \frac{2\pi}{M} k_0 + w_0) ]| \tag{10}
$$

Where the integer DFT index $k_0$ is found to be $k_0 = \text{round}( \frac{M \omega_0}{2\pi} )$. So the peak value of the spectrum is in general not equal to $0.5 A \cdot L$ which you are expecting, which would only be the case when

$$  \frac{2\pi}{M} k_0  =  w_0 \tag{11}$$ for some integer $k_0$ which requires that $\frac{M \omega_0}{2\pi}$ is also an integer, or equivalently

$$ \omega_0 = \frac{ 2\pi }{M}k_0 \tag{12}$$ ,in that special case $|X_L^M[k_0]| = 0.5 A L$ will happen. Otherwise, you will have two complex numbers added and their magnitude is computed.

Note that by selecting $M = L$ when $\omega_0$ meets the criterion (i.e, $ \frac{M \omega_0}{2\pi}$ is an integer), then one can get what is an illusion of lollipop plot indicating only two nonzero samples at the frequencies of $\pm \omega_0$ with magnitude 0.5 $A \cdot L $ and all zeros. But this is a conseqeunce of frequency sampling and not to be taken unnecessarily into too far conclusions.

Here is a MATLAB / OCTAVE code to produce the results

clc; clear all; close all

% Define the sampling grid
Fs = 1000;             % sampling frequency in Hz
Ts = 1/Fs;             % sampling period in seconds

ti = 0;                % initial time (s)
tf = .15;               % final time   (s)
t =  ti :Ts : tf;  % sampling time index t = n*Ts

L = length(t);          % number of samples in x[n] = A*sin(2*pi*f0*n*Ts)
M = L;                  % Number of DFT points to be computed

% Define the sampled signal
A = 0.5;               % amplitude of sine wave

% Signal Frequency Type-1 :
%f0 = 100;              % sine wave frequency in Hz.
%w0 = 2*pi*f0*Ts;        % set an arbitrary frequency
%x = A*cos(w0*[0:L-1]);   % Sine wave sampled at Fs.

% Signal Frequency Type-2:
w0 = 2*pi/M*floor(L/3); % set a special frequency so that DFT PEAK = A/2
x = A*cos(w0*[0:L-1]);   % Sine wave sampled at Fs.

% Calculate DFT of x[n] using FFT function
k0 = ( M*w0/(2*pi));

figure,stem(linspace(-pi,pi,M),fftshift(abs(fft(x,M))/L))
title(['Magnitude of X_L^M[k] for L=',num2str(L),' , M=', num2str(M),...
      ' , A=',num2str(A),' , \omega_0 =', num2str(w0), ' , k_0 = ' num2str(k0)]);
grid on


QPSK constellation from baseband signal My IQ sampler is measuring a QPSK carrier at a set frequency/bandwidth (40MSamp/s), but shifts the signal spectrum to baseband. My signal is 1.5Msym/s, but sampling at 1.5Msym/s gives me a constellation of a circle due to the phase difference between the carrier and the IQ sampler clock.

I still need to synchronize to the carrier, but I now have a 0Hz carrier? Do I try a costas loop with a 0Hz center VCO? So far I haven't had much luck. Any ideas? Thanks everyone.

The following depicts a QPSK Carrier Recovery loop that works on the baseband I and Q samples following proper decisions from a Timing Recovery Loop. For more information on this implementation see

High modulation index PSK - carrier recovery

And for more information on a Timing Recovery implementation that can precede this see:

Symbol timing synchronization using a high sampling rate

(Copying the links here since the question asked would not lead you easily to the prior titles)


Data decomposition method invariant to small shift and scale? Is there a data decomposition method similar to eigenvalue that estimates projection matrix to reduce dimensionality but does not project similar vectors too far away in euclidian distance terms from each other if original data from the same class varies a little in scale, shift and rotation (2D case).

$$y = E  x;$$

e.g. a ECG classification problem example. Cardio cycles have different duration. In addition scale and shift depends on the accuracy on the beat detection. Thus cardio cycles belonging to the same class might be projected to far away due to that variation.

For ECGs, there is work using the Karhunen Lòeve Transform (KLT) (Jager, et al. 1992). I've had success with the KLT technique (Povinelli, 2005).


Relation between Wireless Power Transfer, Wireless Information Transfer and the definition of sum rate I read this paper Charge-then-Forward: Wireless Powered Communication for Multiuser Relay Networks
 ,and it said it on the second page ,and I don't really understand why is the theory ,the author said, right.

  the time of WPT, Wireless Power Transfer, should be as small as possible so that the time for wireless information transmission (WIT) can be maximized for sum-rate maximization

Does anyone know the definition of sum-rate in formula type and the relation of WPT and WIT? I think I can really understand reason if I know these first.

"sum-rate" is just the sum of all rates of communications happening in a network.

The statement is pretty self-evident: The more time your nodes can communicate, the more data they can transfer.


Sinc interpolation in spatial domain I have tried to perform sinc interpolation (in 1D) with the following Matlab code:

Fs=8;
T=1/Fs;
t=0:T:(1-T);
f=1;
x=sin(2*pi*f*t);

up_factor=2;

%% Deduce sinc from Fourier domain

xp=[zeros(1,5) x zeros(1,5)];

Xp=fft(xp);

door1D=abs(xp>0);
sinc1D=fftshift(abs(fft(door1D)));

%plot(door1D);hold on;plot(fftshift(abs(sinc1D)))

%% Interp with sinc in spatial domain
x_up=upsample(x,2);
%plot(x,'b*');hold on;plot(x_up,'r*');
x_up_interp=conv2(x_up,sinc1D,'same');

figure;
plot(x_up_interp./up_factor);
hold on;
plot(x);
hold on;
plot(x_up,'r+');

It seems to work (except for ripples due to the Gibbs phenomenon I guess?). However, I worked this out in a "deductive" manner and I dont completely understand the parameters (period/frequency) of the sinc. The approach was to extract the sinc from the fft of the door function. Then use that sinc as if I had computed it beforehand and convolve the upsampled (not yet interpolated) 1D signal with it.

Can someone help me understand it better? How would I built this sinc ? e.g. from the sinc() function of Matlab. And for it to work, I have convolved with abs of sinc cf. conv2(x_up,sinc1D,'same'); , but this seems strange... can someone explain/develop/correct this ?

REM.:also it is probably badly scaled, but that is another detail

The sinc function actually represents an ideal (brickwall) lowpass filter that's used to complete the interpolation process after the data has been expanded (zero stuffed) properly. So let me outline the time domain approach here:
Assume you have data samples $x[n]$ of length $N$, and you want to upsample this by the integer factor of $L$, yielding a new interpolated data $y[n]$ of length $M = L \times N $ samples.
The first stage is to zero stuff the input $x[n]$; i.e., expand it by $L$ by the expander block :
$$ x[n] \longrightarrow \boxed{ \uparrow L } \longrightarrow x_e[n] $$
where $x_e[n]$ is related to $x[n]$ by the following:
$$ x_e[n] = \begin{cases} { x[n/L] ~~~,~~~ n = 0, \pm L, \pm 2L... \\ ~~~ 0 ~~ ~~~,~~~\text{otherwise} } \end{cases} $$
Then, to complete the interpolation process and fill in the empty (zeroed) samples of $x_e[n]$, one has to lowpass filter $x_e[n]$ by an ideal lowpass filter $h[n]$ with the following frequency domain definition $H(\omega)$ :
$$ H(\omega) = \begin{cases} { ~ L ~ ~~~,~~~ |\omega| < \frac{\pi}{L}  \\ ~ 0 ~ ~~~,~~~\text{otherwise} } \end{cases} $$
The impulse response $h[n]$ of this ideal filter is computed by the inverse discrete-time Fourier transform of $H(\omega)$ and is given by
$$ h[n] = L \frac{ \sin( \frac{\pi}{L} n) } {\pi n} $$
This is an infitely long and non-causal filter, and thus cannot be implemented in this form. (See Hilmar's comments) Practically it's truncated and weighted by a window function, for example by a Hamming or Kaiser window.
The following MATLAB / OCTAVE code represents designing the filter and applying it into data in time domain:
L = 5;      % interpolation factor
N = 500;    % data length

x = hamming(N)'.*randn(1,N);  % generate bandlimited data...

% expanded signal
xe = zeros(1,N*L);
xe(1:L:end) = x;          % generate th expanded signal

% interpolation sinc filter
n = -32:32;               % timing index
h = L*sin(pi*n/L)./(pi*n); % ideal lowpass filter
h(33) = L;                % fill in the zero divison sample
h = hamming(65)'.*h;      % apply weighting window on h[n]

% interpolate:
y = filter(h,1,xe);     % y[n] is the interpolated signal


How to fit Welch's window segments to a time series of variable length I'm trying to write a small Welch's PSD method in R right now, but I've had trouble finding information on how exactly the window segments are fitted to a given time series.

For instance, when the end of the time series does not fall on the end of the last window segment, what are the common ways to deal with this?

I'm assuming the time series is either truncated to perfectly fit the window segments, or zero-padded at the end so another segment can be perfectly fit.  I'm interested in what the standard method is for this.

For example, lets say I have:

A time series of 32 samples
5 window segments, each of 10 samples and 50% overlap

I can either, zero-pad my timeseries until it reaches 35 samples and fit an extra window segment in, or I can truncate my timeseries to 30 samples, which would perfectly fit 5 windows.  Which is the proper course of action here?

I don't think there is one right answer / best practice / standard method in this situation. It would depend on your application. I have seen both of those methods used (zero pad/truncate). Alternatively, you could choose a slightly different window length and/or overlap so that you don't have to truncate too many samples (or pad too many zeroes) near the end of the timeseries.