FFT of a signal using window function

Hi,

I am having the frequency of my signal using the FFT and Welch function. I use the hamming window and rectwin, but the peaks are shifting in both cases as compared to the original signal.I am attaching the original as well as the filtered signal.

Could any one suggest something better?

Kind regards

Tayyaba Bano

[Tayyaba Bano]

I am having the other signal, and highest peak in the signal is getting lower after filetering.

I am not sure and worried whether its filtering out the real frequency value? (please see the following figure) i.e. the value at 0.057 is shifting to a lower peak at 0.05233.

Also, sometimes the filtered value are a little different than the original.i.e. the original 0.33 is not given in filtered , instead at 0.4875(please see the following figure),

Waiting for a kind response

Regards

Tayyaba Bano

[Mathieu NOE]

hello

can you share your data and code ?

[Tayyaba Bano]

Thank you very much for your reply.

Please find the code.

close('all')
Fs = 157;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = length(x);             % Length of signal
t = (0:L-1)*T; 
Y = fft(x-mean(x));
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;

figure()
plot(f,P1)
plot(f,(P1 ./ max(P1)), 'b-', 'LineWidth', 0.35)
xlim([0, 5])
ylim([0, 1.1])
grid on
box on 
xlabel('Frequency (Hz)', 'FontSize',16, FontName='Arial')
ylabel('PSD (dB/Hz)', 'FontSize',16, FontName='Arial')
ax=gca;
ax.XAxis.FontName = 'Cambria Math';
ax.YAxis.FontName = 'Cambria Math';
ax.XAxis.FontSize = 14;
ax.YAxis.FontSize = 14;
ax.GridAlpha = 0.05;
ax.LineWidth = 1.5;

%welches method

noverlap = 0;
nfft=6000;
window= hamming(L);
[pxx,f]= pwelch((x-mean(x)), window, noverlap, nfft, Fs); 
fft_fig = figure();
plot(f,(pxx ./ max(pxx)), 'b-', 'LineWidth', 0.35)
xlim([0, 5])
ylim([0, 1.1])
grid on
box on 

xlabel('Frequency (Hz)', 'FontSize',16, FontName='Arial')
ylabel('PSD (dB/Hz)', 'FontSize',16, FontName='Arial')

Please find original verses filtered signal:

original

filtered:

The value of 0.057 is shifting verses original.

Waiting for any kind response.

Tayyaba

[Mathieu NOE]

hello

in the second case you apply a hamming window on your data

NB that your data is not stationnary , so applying a window may change the amplitude of one peak depending when that frequency appears in the time axis. Should it be strong at the ends and not in the middle of your record , then the window will attenuate it's weight in the fft result

see spectrogram below to visualize how yur signal energy content varies with time

FYI, this is my code , with hanning window (applied on smaller buffer of signal , here 128 samples (after a decimation factor of 10 to focus in the lower frequency range)

clc
clearvars

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename = ('x.xlsx'); 
signal = readmatrix(filename);
signal = detrend(signal);

Fs = 157; % sampling rate 

[samples,channels] = size(signal);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 128;    % 
OVERLAP = 0.99;  % must be between 0 and 0.99

% spectrogram dB scale
spectrogram_dB_scale = 30;  % 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 = 10;
if decim>1
    Fs = Fs/decim;
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),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);
legend(leg_str);
xlim([0 5]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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))) - spectrogram_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)');
    ylim([0 5]);


end

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

[Tayyaba Bano]

Tnakyou very much for your kind help,yes my data is not stationary its changing with time.

So can I not use the filtering with my data?

Regards

Tayyaba Bano

[Mathieu NOE]

why the need to filter ? NB you apply a time window , which is not the same as applying a filter (LP, HP or BP filters)