wlp=0.4999*pi;%下通带截止频率
wls=0.5*pi;%阻带下限频率
wus=0.5*pi;%阻带上限频率
wup=0.5*pi;%上通带截止频率
wc=[(wlp+wls)/2/pi,(wus+wup)/2/pi];%截止频率在通带和阻带边界频率的中点
B=wls-wlp;%迁移域
N=ceil(12*pi/B)-1;%正向取整数,根据过渡带宽等于窗函数主瓣宽求窗函数的最小长度
n=0:N-1;
alpha=(N-1)/2;%求滤波器的相位延迟
m=n-alpha+eps;%eps为MATLAB系统的精度
%hd=sin(wc*m)/(m*pi);%求理想滤波器脉冲响应 %维度不同,无法运行
window=hamming(N);
[h1,w]=freqz(window,1)
figure(1);
subplot(2,1,1)
stem(window,'.');
xlabel('n');
title('海明窗函数');
subplot(2,1,2)
plot(w/pi,20*log(abs(h1)/abs(h1)))
grid on;
xlabel('w/pi');
ylabel('幅度(dB)');
title('海明窗函数频谱');
hn=fir1(N-1,wc,'stop');
[h2,w]=freqz(hn,1,512);
figure(2);
subplot(2,1,1)
stem(n,hn,'.');
xlabel('n');
ylabel('h(n)');
title('海明窗函数单位脉冲响应');
subplot(2,1,2)
plot(w/pi,20*log(abs(h2)/abs(h2(1))));
grid on;
xlabel('w/pi');
ylabel('幅度(dB)');
title('海明窗带阻滤波器的幅度特性');
Ts=0.001;
fs=1/Ts;
x=sin(2*pi*50*n*Ts)+sin(2*pi*125*n*Ts);%原信号
y=x.*window';%加窗函数进行滤波处理
xfft=fft(x,N);
xfft=xfft.*conj(xfft)/N;
y1=fft(y,N);
y2=y1.*conj(y1)/N;
figure(3);%滤除前后的信号对比。
subplot(2,2,1);plot(n,x);grid;
xlabel('时间 (s)');ylabel('幅度');title('输入信号');
subplot(2,2,3);plot(n,y);grid;
xlabel('时间 (s)');ylabel('幅度');title('滤波器输出');
subplot(2,2,2);plot(n*fs/N,xfft);axis([0 fs/2 min(xfft) max(xfft)]);grid;
xlabel('频率 (Hz)');ylabel('Magnitude (dB)');title('输入信号');
subplot(2,2,4);plot(n*fs/N,y2);axis([0 fs/2 min(y2) max(y2)]);grid;
xlabel('频率 (Hz)');ylabel('Magnitude (dB)');title('滤波器输出');