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

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.

Answer 1

[David Goodmanson]

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.

[lekhani gaur]

Here, I have simulated both cases simultaneously. But as the alpha value increases, the spread is less along with the magnitude.

close all;
clear all;
t1 = 0:1e-3:0.5;
w0 = 2*pi*49.5;
alpha = 5;
x1 = 2*cos(w0*t1);
x2 = 2*cos(w0*t1).*exp(-alpha*t1);
plot(t1,x1,'r',t1,x2,'b') 

Ts1 = mean(diff(t1));
fs1 = 1/Ts1;
N1 = length(t1);
n1 = 2^nextpow2(N1);
f1 = fs1*(0:n1-1)/n1;

x1_mags = abs(fft(x1,n1)); 
x2_mags = abs(fft(x2,n1));
 
plot (f1(1:100),x1_mags(1:100),'r',f1(1:100),x2_mags(1:100),'b')

[David Goodmanson]

Hi lekhani,

I am not sure what the point you are trying to make is. I will use y1 (cw) and y2 (exponential decay) here for the freq doman quantities. Certainly y1 looks wider at any given height, but the peak heights are different. If you compare

the width of y1 at 1/sqrt(2) of its height to

the width of y2 at 1/sqrt(2) of its height,

y2 is wider. Since one is coming down by a factor of 1/sqrt(2) from the peak, multiplying the frequency domain function by a constant does not change the width calculation. So another way to do this is to normalize the peaks heights to 1, in which case coming down by 1/sqrt(2) is the same for both and it's easier to see visually. If you plot

x1_magsN = abs(fft(x1,n1))/max(abs(fft(x1,n1))); 
x2_magsN = abs(fft(x2,n1))/max(abs(fft(x2,n1)));
plot (f1(1:100),x1_magsN(1:100),'r',f1(1:100),x2_magsN(1:100),'b')

y2 is wider. [I would have preferred that you not put in the r and b options that reverse the usual order of the plot colors, but there it is]

However, the width of y1 in this case has nothing to do with the parameter alpha. It has to do with an unfavorable time array, which is only made worse by using nextpow2. Take a look at

n = 50000;
delt = .001;
t = (0:n-1)*delt;
delf = 1/(n*delt);
f = (0:n-1)*delf;
f0 = 1;
alpha = .1;
x1 = cos(2*pi*f0*t);
x2 = cos(2*pi*f0*t).*exp(-alpha*t);
figure(1)
plot(t,x1,'r',t,x2,'b')

y1 = fft(x1);
y2 = fft(x2);
figure(2)
% o- option is used to show that y1 is a one-point peak
% multiply y2 by 4 for visual purposes
plot(f,abs(y1),'ro-',f,4*abs(y2),'b')
xlim([0 3])

Here there are exactly 50 cycles in the time window, so y1 has just one point in its peak as it should. y2 is obviously wider.

Compare that to the nextpow2 version:

n = 50000;
delt = .001;
t = (0:n-1)*delt;
m = 2^nextpow2(n)

delf = 1/(m*delt);
f = (0:m-1)*delf;
f0 = 1;
alpha = .1;
x1 = cos(2*pi*f0*t);
x2 = cos(2*pi*f0*t).*exp(-alpha*t);
figure(3)
plot(t,x1,'r',t,x2,'b')

y1 = fft(x1,m);
y2 = fft(x2,m);
figure(4)
% multiply y2 by 4 for visual purposes
plot(f,abs(y1),'r',f,4*abs(y2),'b')
xlim([0 3])

y2 is still wider, but y1 is messed up by the nextpow2 process.

[lekhani gaur]

Thank you so much for such detailed and insightful description. It was really helpful.

[Paul]

Greetings @David Goodmanson,

Do you think that "there is absolutely NO reason" to use nextpow2 in this use case, but there could be a good reason to use it (or zero padding in general) in other use cases?

Or are you recommending to never use zero padding in any analysis?

[David Goodmanson]

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.

Answer 2

[Mathieu NOE]

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

[lekhani gaur]

Thank you for your response. But, If you closely look at the spread of magnitude over frequency range obtained from your code , it is similar to what I got. My question is why is there so much spread of magnitude as highlighted in figure below? Should it not approach to line spectrum with spread over lesser frequency range, say 0.1- 0.2 Hz (ideally no spread should be there).

[Mathieu NOE]

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