implements a digital notch filter to remove the frequency component at 1 Hz.
조회 수: 18 (최근 30일)
이전 댓글 표시
I have a file named X that has data points sampled at 400Hz , time between samples (0.0025sec)
I have been asked to implements a digital notch filter to remove the frequency component at 1 Hz.
using the difference equation:
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)]
where : X(n) my data file X, G=1, p=0.99, w0= 0.0157, n= length(X-1).
I have all the information to implement the equation, but when I use it, it tells me Unrecognized function or variable 'y'.
답변 (1개)
Mathieu NOE
2021년 3월 26일
hello
see example code below
clc
clearvars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% dummy data
Fs = 400;
samples = 2000;
t = (0:samples-1)'*1/Fs;
signal = sin(2*pi*1*t)+0.02*sin(2*pi*50*t)+1e-2*randn(samples,1); % two sine + some noise
%% 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
fc = 1; % notch freq
w0 = 2*pi*fc/Fs;
p = 0.99;
% digital notch (IIR)
num1z=[1 -2*cos(w0) 1];
den1z=[1 -2*p*cos(w0) p^2];
% now let's filter the signal
signal_filtered = filtfilt(num1z,den1z,signal);
figure(1),plot(t,signal,'b',t,signal_filtered,'r');grid
title('Signal');
legend('before notch filter','after notch filter');
xlabel('Time (s)');ylabel(' Amplitude')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 400;
Overlap = 0.75;
[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap); %
sensor_spectrum_dB = 20*log10(fft_spectrum);% convert to dB scale (ref = 1)
% demo findpeaks
[pks,locs]= findpeaks(sensor_spectrum_dB,'SortStr','descend','NPeaks',2);
[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs, NFFT, Overlap);
sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),plot(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')
text(freq(locs)+1,pks,num2str(freq(locs)))
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
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Fourier Analysis and Filtering에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!