1
$\begingroup$

I have hard time in understanding the results I get when I plot the spectrum of the time-derivative of a noise signal:

$|DFT[\frac{dx(t)}{dt}](\omega)|^2$

as a function of frequency $\omega$ (see image below), where $DFT$ is the discrete Fourier Transform, $x(t)$ is a discrete time serie of random samples extracted from a gaussian distribution.

In python:

import numpy as np
import matplotlib.pyplot as plt
x = np.random.randn(100000)
F = np.abs(fft.fft(np.diff(x)))**2
plt.loglog(F[:50000])

The F.transform of a derivative is $i\omega$ times the F.transform. So I'd expect the result to be $\omega^2 |DFT[x]|^2$.

This happens, but what are the features at low frequency? If I run the script several times, a plateau shows up at low frequency, but not always with the same amplitude and not always at the same low frequencies (see in the fig. the differences between various runs of the same script). Probably something strange happens also at high frequency? I heard "leakage" can explain things like this, but I'n not an expert. Can you help me to understand, and correct it?

fft of time derivative

$\endgroup$
2
  • $\begingroup$ plot your data on a linear frequency scale and look at 0 to 40 Hz. The fft is linear in frequency and you’re stretching the low frequency bins on your log frequency plot $\endgroup$
    – user28715
    Commented Oct 30, 2017 at 21:43
  • 1
    $\begingroup$ @StanleyPawlukiewicz Changing scale does not change the result. In loglog, $\omega^2$ is the diagonal line with the right slope in the figure, so loglog is useful. The plateau shows up more clearly in loglog $\endgroup$
    – scrx2
    Commented Oct 31, 2017 at 8:16

2 Answers 2

2
$\begingroup$

I think your python code is wrong. Below I put a modified and corrected MATLAb/OCTAVE code instead. I hope it's clear from the code about the difference.

 N = 32*256;                       % Length of the time domain signal
 x = randn(1,N);                   % Normal type random vector
 wk = linspace(-pi,pi-2*pi/N,N);   % DFT frequencies from -pi to pi
 X = fftshift(fft(x,N));           % DFT X[k] of x[n] FROM -pi to pi
 Px = (wk.^2).*(abs(X).^2);        % Simulated |DFT|^2 of the derivative signal  
 figure,plot(linspace(-1,1,N),Px)  % plot the $DFT|^2 from -pi to pi
 title('Simplified Periodogram');  % QUADRATIC behaviour is recognised.
 xlabel('Normalized Frequency (x \pi rad/sample)'); % No DC plateau ?

 % LOG-LOG plot.
 figure,loglog(10E3*wk(N/2+1:end),Px(N/2+1:end))

 % Alternatively: compute PSD based on time-domain signal
 pX = fftshift(j*wk.*X);            % DFT of the derivative
 dx = real(ifft(pX,N));                % time-domain signal from I-DFT
 figure,periodogram(dx,hamming(N));    % compute & display its Periodogram
 figure,pwelch(dx);                    % compute & display its Welch

The resulting plots are:

The loglog plot of MATLAB is the following:

enter image description here

there is no DC plateau here. I believe it's your diff function which adds the DC content to your differentiated data. You should apply the DSP standard diagnostic procedure to investigate your case and find out which step deviates from the expected behaviour...

$\endgroup$
6
  • $\begingroup$ but it seems you take the fft(derivative) already as $\omega^2 fft(x)^2$. I would like to know why (fft(diff(x)))^2 is not equal to $\omega^2 fft(x)^2$ at low frequency $\endgroup$
    – scrx2
    Commented Oct 31, 2017 at 8:20
  • $\begingroup$ They should be equal but where is (fft(diff(x)))^2 inyour code ? $\endgroup$
    – Fat32
    Commented Oct 31, 2017 at 9:20
  • $\begingroup$ ouch, sorry! the code is edited $\endgroup$
    – scrx2
    Commented Oct 31, 2017 at 13:36
  • $\begingroup$ Hmm now you have a different code... So the results must be the same provided that the function diff(x) indeed returns the samples of the differentiated time domain signal. The only mismatch is due to the display plot funciton being loglog. (I assume that your random number generator is also compatbile with Matlab's...) $\endgroup$
    – Fat32
    Commented Oct 31, 2017 at 15:38
  • 1
    $\begingroup$ if it's not the loglog function then it could be the diff function. As there's nothing else left for a malfunction for this simple code. $\endgroup$
    – Fat32
    Commented Oct 31, 2017 at 22:38
0
$\begingroup$

Try applying a window function (e.g. Hann). It really seems that the fft leakage of a violet noise signal using a rectangular (or no) window somehow results in a flat plateau (white noise) below some frequency. I also thought about aliasing effects, but the Hann window wouldn't remove aliases, right?

Do you know ALANOISE? If you set it to "white PM", 40000 points, and execute "Generate" and "Compute", you will see a plateau, which is roughly 7 y-units above the apparent PSD value at x-value 0. If you repeat "Generate" and "Compute", the plateau level varies quite a lot (I had instances with no visible plateau), but on average it is about 7 units up. If you save that signal and import it to Matlab/Octave/Scilab, a (rectangular windowed) fft will show the same plateau, while a Hann windowed fft will show a (relatively) continuous slope. I'd appreciate if someone could give an explanation on this effect.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.