利用频率采取法设计50HZ陷波器（matlab代码）

%% FIR陷波器设计
%% -----------------------频率采取法Frequency adoption Method--------------------%%
% 目的：设计一个陷波器阻带在50±1.5Hz以内,采样频率为400Hz的滤波器，
% 并要求通带最大衰减为0.1dB，阻带最小衰减为80dB。
clc;
clear;close all;
rp=0.1; % 通带波纹最大衰减为0.1dB
rs=80; % 阻带衰减为80dB
fs=400; % 采样频率为400Hz
Ts=1/fs;
% fp1=48.5;fp2=51.5;
% fs1=49.9;fs2=50.1;
wp_l=48.52.pi/fs;wp_h=51.52.pi/fs; % 根据阻带：49.9Hz_{50.1Hz，通带：0}48.5Hz和51.5Hz~pi，w=mo_wT=mo_w/fs=2pif/fs
ws_l=49.62.pi/fs;ws_h=50.42.*pi/fs; % 确定两个通带截止频率和两个阻带截止频率。
Bt=ws_l-wp_l; % 计算过渡带宽度
m=2; % 过渡点个数m=1
N=ceil((m+1)*2.*pi/Bt)+1; % 估计采样点N
alph=(N-1)/2;
k0=(0:N-1);
wk=2.*pi/N;
T1=0.11; T2=0.59;% 过度采样样本点
hk=[ones(1,ceil(wp_l/wk)),T2,T1,zeros(1,ceil((Bt/wk)-1)),T1,T2,ones(1,ceil((1.4775.pi/wk)+5)),T2,T1,zeros(1,ceil((Bt/wk)-1)),T1,T2,ones(1,ceil(wp_l/wk))]; % 幅度采样向量（偶对称）
wdl=[0 wp_l/pi ws_l/pi ws_h/pi wp_h/pi 1];
k1=0:floor((N-1)/2); k2=floor((N-1)/2)+1:N-1;
angh=[-alph(2.*pi)k1/N,alph(2.pi)(N-k2)/N]; % 相位采样向量
hdk=hk.exp(jangh);% 构成频率采样向量hdk
hn=real(ifft(hdk));% hn=IDFT[hk], real只取实部，忽略计算误差引起的虚部
[hw,w]=freqz(hn,1); %求hn的频谱

%% =波形绘制=%%
figure;
% stem(k0.wkfs/(2pi)/4,hn,’.’);
% title(‘impulse response’); xlabel(‘n’);ylabel(‘h(n)’);grid;
% axis([0,max(k0.wkfs/(2pi)/4),0,max(hn)]);
% grid on;
subplot(3,2,1);
plot(k0.wk.fs/(2pi),abs(hdk),’.b-’);
title(‘amplitude response’); xlabel(‘f/Hz’);ylabel(‘H(e^jf)’);grid;
axis([40,60,0,1.1]); % 定义横坐标和纵坐标的范围
grid on;
subplot(3,2,2);
plot(wfs/(2pi),20log10(abs(hw))/max(abs(hw)));
title(‘损耗函数’); xlabel(‘f/Hz’);ylabel(‘A(f)/dB’);
axis([45,55,min(20log10(abs(hw))/max(abs(hw))),10]); % 定义横坐标和纵坐标的范围
grid on;
%% =信号产生=%%
Nlen=2038400;
n=1:Nlen-1;
xn=sin(2pi48nTs)+sin(2pi48.5nTs)+sin(2pi49nTs)+sin(2pi50.nTs)+sin(2pi51nTs)+sin(2pi51.5nTs)+sin(2pi52nTs); % 信号叠加
k=2(0:Nlen-1)/Nlen;
xnk=Tsfft(xn,Nlen);
yn=fftfilt(hn,xn,Nlen);
yk=fft(yn,Nlen); % 计算平率响应函数：DFT[hn]yn_boxcar=fftfilt(hn_1,xn,Nlen);
%%=波形绘制=%%
subplot(3,2,3);
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(‘幅值’);
axis([0,1000,min(xn),max(xn)]); % 定义横坐标和纵坐标的范围
subplot(3,2,4);
plot(kfs/2,abs(xnk),‘b’); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2
title(‘数字滤波后输出信号x(n)频谱特性’); xlabel(‘f/Hz’);ylabel(‘H(e^jf)’);
axis([45,55,0,max(abs(xnk))]); % 定义横坐标和纵坐标的范围
subplot(3,2,5);
plot(n,yn); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=kfs/2
title(‘滤波后输出信号y(n)波形图’); xlabel(‘f/Hz’);ylabel(‘幅值’);
axis([0,1000,min(yn),max(yn)]); % 定义横坐标和纵坐标的范围
subplot(3,2,6);
plot(kfs/2,abs(yk),‘b’); % k=w/pi,w=mo_wT=mo_w/fs,mo_w=2pif;得到f=k*fs/2
title(‘频率采样数字滤波后输出信号y(n)频谱特性’); xlabel(‘f/Hz’);ylabel(‘H(e^jf)’);
axis([45,55,0,max(abs(yk))]); % 定义横坐标和纵坐标的范围