matlab陷波器函数,利用窗函数设计陷波器（matlab代码）

%% FIR陷波器设计

%% =窗函数设计法window function Method========%%

% 目的：设计一个陷波器阻带在50±1.5Hz以内,采样频率为400Hz的滤波器，

% 并要求通带最大衰减为0.1dB，阻带最小衰减为80dB。

clc;

clear;close all;

rp=0.1; % 通带波纹最大衰减为0.1dB

rs=80; % 阻带衰减为80dB

fs=400; % 采样频率为400Hz

wp_l=48.52pi/fs;wp_h=51.52pi/fs; % 根据阻带：49.9Hz50.1Hz，通带：048.5Hz和51.5Hz~pi，w=mo_wT=mo_w/fs=2pif/fs

ws_l=49.92pi/fs;ws_h=50.12pi/fs; % 确定两个通带截止频率和两个阻带截止频率。

Bt=ws_l-wp_l; % 计算过渡带宽度

wc=[(wp_l+ws_l)/(2pi),(wp_h+ws_h)/(2*pi)];% 截止频率在通带和阻带边界频率的中点;(关于π归一化)

%% ===窗函数设计=%%

% 1.1 矩形窗

N_1_0=ceil(1.8*pi/Bt);

N_1=N_1_0+mod(N_1_0+1,2)

hn_1=fir1(N_1-1,wc,‘stop’,boxcar(N_1));

[hw_1,w_1]=freqz(hn_1);% 求hn的频谱

% 1.2 三角窗

N_2_0=ceil(6.1*pi/Bt);

N_2=N_2_0+mod(N_2_0+1,2)

hn_2=fir1(N_2-1,wc,‘stop’,bartlett(N_2));

[hw_2,w_2]=freqz(hn_2);% 求hn的频谱

% 1.3 汉宁窗

N_3_0=ceil(6.2*pi/Bt);

N_3=N_3_0+mod(N_3_0+1,2)

hn_3=fir1(N_3-1,wc,‘stop’,hanning(N_3));

[hw_3,w_3]=freqz(hn_3);% 求hn的频谱

% 1.4 哈明窗

N_4_0=ceil(6.6*pi/Bt);

N_4=N_4_0+mod(N_4_0+1,2)

hn_4=fir1(N_4-1,wc,‘stop’,hamming(N_4));

[hw_4,w_4]=freqz(hn_4);%求hn的频谱

% 1.5 布莱克曼窗

N_5_0=ceil(11*pi/Bt);

N_5=N_5_0+mod(N_5_0+1,2)

hn_5=fir1(N_5-1,wc,‘stop’,blackman(N_5));

[hw_5,w_5]=freqz(hn_5); % 求hn的频谱

% 1.6 凯塞窗(=7.865)

alph=0.5842*(rs-21)^0.4+0.07886*(rs-21);% 根据(7.2.16)式计算Kaiser窗的控制参数a

N_6_0=ceil((rs-8)/(2.286*Bt)); % 根据(7.2.17)式计算Kaiser窗的所需阶数N

N_6 = N_6_0+mod(N_6_0+1,2)

hn_6=fir1(N_6-1,wc,‘stop’,kaiser(N_6,alph));

[hw_6,w_6]=freqz(hn_6); % 求hn的频谱

%% =各个滤波器的幅频响应对比分析=%%

A1 = w_1fs/(2pi);

A2 = w_2fs/(2pi);

A3 = w_3fs/(2pi);

A4 = w_4fs/(2pi);

A5 = w_5fs/(2pi);

A6 = w_6fs/(2pi);

figure; % 画图

subplot(2,1,1);

plot(A1,abs(hw_1),A2,abs(hw_2),A3,abs(hw_3),A4,abs(hw_4),A5,abs(hw_5),A6,abs(hw_6));grid;

title(‘频率响应幅度’); xlabel(‘f／Hz’); ylabel(‘H(ejf)’);

legend(‘boxcar’,‘bartlett’,‘hanning’,‘hamming’,‘blackman’,‘Kaiser’);

axis([47,53,0,1.5]);% 定义横坐标和纵坐标的范围

subplot(2,1,2);

plot(A1,20log10(abs(hw_1))/max(abs(hw_1)),A2,20log10(abs(hw_2))/max(abs(hw_2)),A3,20log10(abs(hw_3))/max(abs(hw_3)),A4,20log10(abs(hw_4))/max(abs(hw_4)),A5,20log10(abs(hw_5))/max(abs(hw_5)),A6,20log10(abs(hw_6))/max(abs(hw_6)));

title(‘损耗函数’); xlabel(‘f/Hz’);ylabel(‘A(f)/dB’);

legend(‘boxcar’,‘bartlett’,‘hanning’,‘hamming’,‘blackman’,‘Kaiser’);

axis([47,53,min(20*log10(abs(hw_6))/max(abs(hw_6))),10]); % 定义横坐标和纵坐标的范围

grid on;

%% =信号产生=%%

Nlen=2038400;

n=1:Nlen-1;

Ts=1/fs;

xn=sin(2pi48nTs)+sin(2pi48.5nTs)+sin(2pi49nTs)+sin(2pi50.nTs)+sin(2pi51nTs)+sin(2pi51.5nTs)+sin(2pi52nTs); % 信号叠加

k=2*(0:Nlen-1)/Nlen;

yn_boxcar=fftfilt(hn_1,xn,Nlen); % 复合信号通过boxcar滤波器后的信号yn;

yn_bartlett=fftfilt(hn_2,xn,Nlen); % 复合信号通过bartlett滤波器后的信号yn;

yn_hanning=fftfilt(hn_3,xn,Nlen); % 复合信号通过hanning滤波器后的信号yn;

yn_hamming=fftfilt(hn_4,xn,Nlen); % 复合信号通过hamming滤波器后的信号yn;

yn_blackman=fftfilt(hn_5,xn,Nlen); % 复合信号通过blackman滤波器后的信号yn;

yn_kaiser=fftfilt(hn_6,xn,Nlen); % 复合信号通过kaiser滤波器后的信号yn;

xnk_T_fft=Tsfft(xn,Nlen); % 对信号做len点FFT变换

yk_boxcar=Tsfft(yn_boxcar,Nlen); % 对信号做M点FFT变换

yk_bartlett=Tsfft(yn_bartlett,Nlen); % 对信号做M点FFT变换

yk_hanning=Tsfft(yn_hanning,Nlen); % 对信号做M点FFT变换

yk_hamming=Tsfft(yn_hamming,Nlen); % 对信号做M点FFT变换

yk_blackman=Tsfft(yn_blackman,Nlen); % 对信号做M点FFT变换

yk_kaiser=Tsfft(yn_kaiser,Nlen); % 对信号做M点FFT变换

%% =波形绘制=%%

figure; % 画图

subplot(4,2,1);

plot(n,xn); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2

title(‘输入信号x(n)原波形图’); xlabel(‘f/Hz’);ylabel(‘H(ejf)’);

axis([0,4000,min(xn),max(xn)]); % 定义横坐标和纵坐标的范围

subplot(4,2,2);

plot(kfs/2,abs(xnk_T_fft),‘b’); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2

title(‘数字滤波后输出信号x(n)频谱特性’); xlabel(‘f/Hz’);ylabel(‘H(ejf)’);

axis([45,55,0,max(abs(xnk_T_fft))]); % 定义横坐标和纵坐标的范围

subplot(4,2,3);

plot(kfs/2,abs(yk_boxcar),‘b’); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2

title(‘boxcar滤波后输出信号y(n)频谱特性’); xlabel(‘f/Hz’);ylabel(‘H(ejf)’);

axis([45,55,0,max(abs(yk_boxcar))]); % 定义横坐标和纵坐标的范围

subplot(4,2,4);

plot(kfs/2,abs(yk_bartlett),‘b’); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2

title(‘bartlett滤波后输出信号y(n)频谱特性’); xlabel(‘f/Hz’);ylabel(‘H(ejf)’);

axis([45,55,0,max(abs(yk_bartlett))]); % 定义横坐标和纵坐标的范围

subplot(4,2,5);

plot(kfs/2,abs(yk_hanning),‘b’); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2

title(‘hanning滤波后输出信号y(n)频谱特性’); xlabel(‘f/Hz’);ylabel(‘H(ejf)’);

axis([45,55,0,max(abs(yk_hanning))]); % 定义横坐标和纵坐标的范围

subplot(4,2,6);

plot(kfs/2,abs(yk_hamming),‘b’); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2

title(‘hamming滤波后输出信号y(n)频谱特性’); xlabel(‘f/Hz’);ylabel(‘H(ejf)’);

axis([45,55,0,max(abs(yk_hamming))]); % 定义横坐标和纵坐标的范围

subplot(4,2,7);

plot(kfs/2,abs(yk_blackman),‘b’); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2

title(‘blackman滤波后输出信号y(n)频谱特性’); xlabel(‘f/Hz’);ylabel(‘H(ejf)’);

axis([45,55,0,max(abs(yk_blackman))]); % 定义横坐标和纵坐标的范围

subplot(4,2,8);

plot(kfs/2,abs(yk_kaiser),‘b’); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2

title(‘kaiser滤波后输出信号y(n)频谱特性’); xlabel(‘f/Hz’);ylabel(‘H(ejf)’);

axis([45,55,0,max(abs(yk_kaiser))]); % 定义横坐标和纵坐标的范围