Why is the magnitude spectrum spreaded over certain frequency range instead of line spectrum in FFT?

18 views (last 30 days)
I have obtained a displacement time history of a system using ODE45 solver. I need to get the frequency content from this time domain data. It is apparent after looking at the time domain data that there should be two frequencies. However, the magnitude spectrum obtained from FFT shows spread of magnitude over a frequency range of ~ 4 Hz rather than line spectrum at two frequencies. I have attached the time doain data in excel file and . mlx file for the MATLAB code. The ODE45 solver yields variable time step, hence I have resampled the data before doing FFT analysis. Also, I tried FFT with time doain data obtained by fixing the time step in ODE45. But in both cases FFT response are similar.
What is the posible reason of it? I have also checked with windowing and zero padding. Zero padding only increses the frequency resolution, but the spread of magnitude over frequency range of ~ 4 Hz is still there.
Kindly suggest if I am missing anything here.
Thank you.

Accepted Answer

David Goodmanson
David Goodmanson on 14 Jan 2022
Edited: David Goodmanson on 14 Jan 2022
Hi lekhani,
The spread is happening because you don't have a pure, CW, constant-amplitude cosine wave in the time domain. That's what it would take to get a sharp peak. Rather you have a wave that is decaying, whose amplitude is decreasing roughly as e^(-alpha*t). There are two frequencies involved but for simplicity let's take just the higher, 49.5 Hz one. An additional second frequency does not change the argument.
A reasonably good one-frequency analytic fit for the signal is
w0 = 2*pi*49.5;
alpha = 2.5
xana = .035*cos(w0*t1).*exp(-alpha*t1)
as you can see by eye with
figure(1); grid on
w0 = 2*pi*49.5;
alpha = 2.5
xana = .035*cos(w0*t1).*exp(-alpha*t1)
plot(t1,x1,t1,xana)
However, after your fft, the frequency domain peak appears at about 38.5 Hz. This is a SERIOUS error
and it arises from using
n1 = 2^nextpow2(N1)
f1 = fs1*(0:n1-1)/n1;
That is, calculating a new frequency array based on a different number of points than was used for the fft. Now the frequency array and the time array have a different number of points, which for the fft should not be. This is a variation on the countless number of ways that nextpow2 can cause harm.. I will say a bit more about nextpow2 at the end but I strongly suggest that you DON'T use it. And tell your friends. For now if you replace
n1 = 2^nextpow2(N1)
with
n1 = N1
and the same for N2, the frequencies appear where they should be. No need for powers of 2.
The fourier transform of xana, which I will call yana, has peaks at positive and negative frequencies. The positive frequency peak is to a very good approximation
yana = 1/(i(w-w0)+alpha)
with magnitude
abs(yana) = 1/sqrt( (w-w0)^2+alpha^2 )
which is not a super-sharp peak about w0, but one with a nonzero width that depends on the value of alpha.
The standard way to calculate the width is as follows.
When w = w0 the max height of the peak is 1/alpha. If you go to the nearby frequencies on each side, w = w0 +- alpha, the value of abs(yana) is down from the peak height by a factor of 1/sqrt(2) as you can check. At those locations the full width of the peak is (w0+alpha)-(w0-alpha) = 2*alpha. That's the width in the w domain, circular frequency. Since f = w/(2*pi), the width in the regular frequency domain is 2*alpha/(2*pi) = alpha/pi.
So there is a prediction about the real fft data: dropping down from the peak by a factor of 1/sqrt(2), the full peak width should be alpha/pi. The data itself has sets x1 and x2, with x2 being almost the same as x1 but with a shorter time record. So I will concentrate on x1_mags, the blue curve in figure 2. The height is 37.7, and after dropping down by sqrt(2) to 27.7, the width at that height should be 2.5/pi = 0.8 Hz. If you zoom in on a the height 27.7 the full width is actually 0.7 Hz. Fairly close.
As far as nextpow2, the usual process would be
N1 = length(t1)
n1 = 2^nextpow2(N1)
y1 = fft(x1,n1) % pads x1 with zeros to length n1
f1 = fs1*(0:n1-1)/n1; % as you did
but there is absolutely NO reason to do this. There is supposedly a gain in speed with use of 2^(some power) but in fact there is no appreciable gain in speed (except for exceptionable circumstances), padding the signal generally messes up the frequency spectrum and there are many opportunities for error.
  5 Comments
David Goodmanson
David Goodmanson on 20 Jan 2022
Hello Paul,
There appear to be at least a couple of favorable situations for zeropadding, but I cannot think of a reason to use nextpow2 for that or anything else.
One use of zeropadding is that it provides a finer array with smaller delta_f in the frequency domain. This is most often done by increasing the number of time points by some integer factor. Then the original frequencies are still part of the new frequency array and comparison is easy. That definitely does not happen with nextpow2. The new array has strange spacing compared to the old one, so using nextpow2 is a disincentive. (One exception is if the original array length is a power of 2; then increasing the array length some power of 2 works. You wouldn't do that using nextpow2, though).
The speed of an fft basically boils down to the prime factors of the array length N. if [1] N is a large prime number, or even if [2] N has some medium size prime factors, the fft is slow. A speed advantage occurs when N is the product of small primes, even if there are a lot of them. For N as in [2], zeropadding to 2^n provides a significant speed advantage, typically on the order of a factor of 10. But 2^n is not the only case around.
Arrays with N = 1,2,3,4,5,7*!0^m (six possibilities) also have small prime factors, and on my pc with N in the range from 1e3 to 3e8 the fft actually runs faster than it does when those lengths are zeropadded to the next 2^n. In some cases, constant*10^m even runs faster than the next smaller case of 2^n. So I see the supposed supremacy of 2^n as basically, if not a myth, then unjustifiable folklore that 'everybody knows about'.
This result might be hardware-dependent, but I am using your basic windows 10 intel cpu so I think many people would get similar results.
Zeropadding will always have consequences, but zeropadding to 1,2,4,5*10^m has another advantage over 2^n than just speed. Suppose the time array has a convenient array spacing such as delta_t = 1e-6. When the array is zeropadded to length N, then by the golden rule for the fft,
delta_t*delta_f = 1/N.
If N = 1,2,4,5*10^m then delta_f will be a fairly convenient quantity. But if you zeropad, say 7800 points to 2^n = 8192, the resulting delta_f is the rather meaningless 122.0703125 Hz. I am not saying that this is a major advantage of 1,2,4,5*10^m over 2^n, but it doesn't hurt.

Sign in to comment.

More Answers (1)

Mathieu NOE
Mathieu NOE on 13 Jan 2022
hello
I took the 3rd data and used my usual tool and could get the two frequencies
clc
clearvars
option_notch = 0; % 0 = without notch filter , 1 = with notch filter
fc_notch = 50; % notch freq
option_LPF = 0; % 0 = without low pass filter , 1 = with low pass filter
fc_lpf = 50; % LPF cut off freq
N_lpf = 2; % filter order
option_HPF = 0; % 0 = without high pass filter , 1 = with high pass filter
fc_hpf = 150; % HPF cut off freq
N_hpf = 2; % filter order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
out = importdata('data3.txt');
data = out.data;
time = data(:,1);
signal = data(:,2);
dt = mean(diff(time));
Fs = 1/dt;
[samples,channels] = size(signal);
% recreate time vector for perfect time accuracy
time = (0:samples-1)*dt;
%% notch filter section %%%%%%
% y(n)=G*[x(n)-2*cos(w0)*x(n-1)+x(n-2)]+[2*p cos(w0)*y(n-1)-p^2 y(n-2)]
% this difference equation can be converted to IIR filter numerator /
% denominator
if option_notch ~= 0
w0 = 2*pi*fc_notch/Fs;
p = 0.995;
% digital notch (IIR)
num1z=[1 -2*cos(w0) 1];
den1z=[1 -2*p*cos(w0) p^2];
% now let's filter the signal
signal = filtfilt(num1z,den1z,signal);
end
%% low pass filter section %%%%%%
if option_LPF ~= 0
w0_lpf = 2*fc_lpf/Fs;
% digital notch (IIR)
[b,a] = butter(N_lpf,w0_lpf);
% now let's filter the signal
signal = filtfilt(b,a,signal);
end
%% high pass filter section %%%%%%
if option_HPF ~= 0
w0_hpf = 2*fc_hpf/Fs;
% digital notch (IIR)
[b,a] = butter(N_hpf,w0_hpf,'high');
% now let's filter the signal
signal = filtfilt(b,a,signal);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 5000; %
OVERLAP = 0.95;
% dB scale
dB_scale = 120; % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 4;
if decim>1
for ck = 1:channels
newsignal(:,ck) = decimate(signal(:,ck),decim);
Fs = Fs/decim;
end
signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;
%%%%%% legend structure %%%%%%%%
for ck = 1:channels
leg_str{ck} = ['Channel ' num2str(ck) ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),plot(time,signal);grid on
title(['Time plot / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend(leg_str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);
% apply A weigthing if needed
if option_w == 1
pondA_dB = pondA_function(freq);
sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
my_ylabel = ('Amplitude (dB (A))');
else
my_ylabel = ('Amplitude (dB (L))');
end
figure(2),plot(freq,sensor_spectrum_dB);grid on
df = freq(2)-freq(1); % frequency resolution
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
mm = 5*ceil(max(sensor_spectrum_dB,[],'all')/5); % 5 dB rounded max level
ylim([mm-dB_scale mm]);
legend(leg_str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
[sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg)); % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% apply A weigthing if needed
if option_w == 1
pondA_dB = pondA_function(fsg);
sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
my_title = ('Spectrogram (dB (A))');
else
my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range :
% saturation_dB = 60; % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2+ck);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid on
df = fsg(2)-fsg(1); % freq resolution
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
xlabel('Time (s)');ylabel('Frequency (Hz)');
end
sound(signal,Fs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pondA_dB = pondA_function(f)
% dB (A) weighting curve
n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
pondA = n./r;
pondA_dB = 20*log10(pondA(:));
end
function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).
% Linear averaging
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,channels);
s_tmp((1:samples),:) = signal;
signal = s_tmp;
samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
% compute fft with overlap
offset = fix((1-Overlap)*nfft);
spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
% % for info is equivalent to :
% noverlap = Overlap*nfft;
% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows
% main loop
fft_spectrum = 0;
for i=1:spectnum
start = (i-1)*offset;
sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only
end
fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum % Select first half
if rem(nfft,2) % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end
  2 Comments
Mathieu NOE
Mathieu NOE on 14 Jan 2022
hi
in order to have "no spread" in the fft result, the frequencies of your signal must be matching exactly the frequency bins of the fft (which is discrete) , which happens rarely as most of the time the signal frequencies are unknown

Sign in to comment.

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!