How to find Frequencies to filter?

Prashant Kumar answered . 2024-12-21 20:02:02

clc
clear all

% hi 
% may I Jump in the conversation ? 
% there are so many questions in this forum dealing with fft , and many times it is not always clear for everybody how to get the best result according to what kind of signals are being analysed (random, harmonics, stationnary or not ...).
% I like to propose this little demo based on pwelch and spegram to give another point of view and also let the people play with fft size, overlap , windows and see the effect on frequency resolution, amplitude resolution etc...


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dummy signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fs = 5000;
samples = 255001;
dt = 1/Fs;

t = (0:dt:(samples-1)*dt);
omega = 2*pi*(25+20*t);

sensordata = randn(size(t)) + 5*sin(omega.*t);  % signal = linear chirp + noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% FFT parameters
samples = length(sensordata);
NFFT = 512;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% option 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sensor_spectrum, freq] = pwelch(sensordata,w,NOVERLAP,NFFT,Fs);

figure(1),plot(freq,20*log10(sensor_spectrum));
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude (dB)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% option 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

saturation_dB = 80;  % dB range scale (means , the lowest displayed level is 80 dB below the max level)
% fmin = 1;
% fmax = Fs/2;

[sg,fsg,tsg] = specgram(sensordata,NFFT,Fs,w,NOVERLAP);  

% 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

% saturation to given dB range
mini_dB = round(max(max(sg_dBpeak))) - saturation_dB;
sg_dBpeak(sg_dBpeak<mini_dB) = mini_dB;

% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);axis('xy');colorbar('vert');
title(['Spectrogram  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');

Not satisfied with the answer ?? ASK NOW