% Design a bandstop filter
% see doc fir1
% b=fir1(n,Wn) returns row vector b containing the n+1 coefficients of an order n lowpass FIR filter.
% This is a Hamming-window based, linear-phase filter with normalized cutoff frequency Wn.
% Wn is a number between 0 and 1, where 1 corresponds to the Nyquist frequency.
% b = fir1(n,Wn,'ftype') specifies a filter type, where 'ftype' is:
% 'stop' for a bandstop filter, if Wn=[w1 w2]. The stopband frequency range is specified by this interval.
clc
clear all
n = 40; % filter length
T = 0.04; % sampling time
sampfreq=1/T; % sampling frequency
wc1 = 8*pi; % cutoff frequency
wc2 = 12*pi; % cutoff frequency
Omegac1 = wc1/(0.5*sampfreq*2*pi); % normalized cutoff frequency
Omegac2 = wc2/(0.5*sampfreq*2*pi); % normalized cutoff frequency
Omegac = [Omegac1,Omegac2];
b = fir1(n, Omegac,'stop'); % create a FIR digital filter
[hd, tn] = impz(b, 1, 150); % impulse response
t = [0:150]*T;
x = 1 + 2*cos(6*pi*t) + 3*cos(10*pi*t); % 滤波前的信号
y = filter(b, 1, x); % 对 x 进行滤波
sigLen = length(x);
NFFT = 2^nextpow2(sigLen); % Next power of 2 from length of x
X = fft(x,NFFT)/sigLen;
Y = fft(y,NFFT)/sigLen;
f = sampfreq/2*linspace(0,1,NFFT/2+1);
figure
set(gcf,'color','w');
subplot(221)
plot(t,x)
xlabel('time(s)');
ylabel('original signal');
subplot(222)
plot(f,2*abs(X(1:NFFT/2+1)))
ylabel('|X(f)|');
xlabel('frequency(Hz)');
subplot(223)
plot(t,y)
ylabel('filtered signal');
xlabel('time(s)');
subplot(224)
plot(f,2*abs(Y(1:NFFT/2+1)))
ylabel('|Y(f)|');
xlabel('frequency(Hz)');
