NFFT = 1024*4;
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);
Fs = 1000;
t = 0:1/Fs:10-1/Fs;
signal = cos(2*pi*50*t)+cos(2*pi*100*t);
fc = 50;
wc = 2*pi*fc;
Q = 2;
w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);
b0 = 1;
b1 = -2*cos(w0);
b2 = 1;
a0 = 1 + alpha;
a1 = -2*cos(w0);
a2 = 1 - alpha;
num1=[1/wc^2 0 1];
den1=[1/wc^2 1/(wc*Q) 1];
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
freq = linspace(fc-1,fc+1,200);
[g1,p1] = bode(num1,den1,2*pi*freq);
g1db = 20*log10(g1);
[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1);
figure(1);
plot(freq,g1db,'b',freq,gd1db,'+r');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
legend('analog','digital');
xlabel('Fréquence (Hz)');
ylabel(' dB')
signal_filtered = filtfilt(num1z,den1z,signal);
[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs);
sensor_spectrum_dB = 20*log10(sensor_spectrum);
[sensor_spectrum_filtered, freq] = pwelch(signal_filtered,w,NOVERLAP,NFFT,Fs);
sensor_spectrum_filtered_dB = 20*log10(sensor_spectrum_filtered);
figure(2),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')
Best Answer