帮忙设计一个基于频率抽样法的FIR高通滤波器的设计与仿真 用MATLAB，要具体程序和仿真步骤，参数 非常感谢

用频率采样设计法设计一个线性相位FIR高通滤波器，其理想特性为 

close all;clear all;clc;

N=33;wc=pi/3;

N1=fix(wc/(2*pi/N));N2=N-2*N1-1;

A=[zeros(1,N1+1),ones(1,N2),zeros(1,N1)]; 

A(N1+2)=0.3904;A(N-N1)=0.3904;

theta=-pi*[0:N-1]*(N-1)/N;

H=A.*exp(j*theta); 

h=real(ifft(H));v=1:N;

subplot(2,2,1),

plot(v ,A,'k*');title('频率样本');ylabel('H(k)');axis([0,fix(N*1.1),-0.1,1.1]);

subplot(2,2,2),

stem(v ,h,'k');title('脉冲响应');ylabel('h(n)');

axis([0,fix(N*1.1),min(h)*1.1,max(h)*1.1]);

M=500;

nx=[1:N];

w=linspace(0,pi,M);

X=h*exp(-j*nx'*w);

subplot(2,2,3),

plot(w./pi,abs(X),'k');xlabel('\omega/\pi');ylabel('Hd(w)');

title('幅度响应');axis([0,1,-0.1,1.3]);

subplot(2,2,4),

plot(w./pi,20*log10(abs(X)),'k');title('幅度响应'),

xlabel('\omega/\pi');ylabel('dB');axis([0,1,-80,10]);

clear all;close all;clc;
N=51; %滤波器系数的长度（时域点数,可以自定义，一般为奇数）
M=N-1; %辅助参数
Wc=pi/4; %理想截止频率
%符号运算，求hn→FIR-DLP的时域序列
syms k n real %指定符号变量→w:频率，n：离散时间
%N为偶数时
Hd_Even_Front=exp(-j*k*pi*(N-1)/N); %0,1,...,N/2-1
Hd_Even_Half=0; %N/2
Hd_Even_Back=-exp(-j*k*pi*(N-1)/N); %N/2+1,N/2+2,...,N-1
%N为奇数时
Hd_Odd=exp(-j*k*pi*(N-1)/N); %0,1,...,N-1
%-----------------------------------------------------------------
%先计算Hd,再计算hn→IDFT(Hd)
nc=floor(N*Wc/(2*pi)); %计算通带Hd截止(Matlab)坐标
Hd_value=0.3830; %改进滤波器性能参数(0.3～0.5)
if 0~=rem(N,2) %N为奇数
for i=0:N-1
if i<=nc
Hd(i+1)=double(limit(Hd_Odd,k,i)); %计算Hd_奇数
elseif i>=N-nc
Hd(i+1)=double(limit(Hd_Odd,k,i)); %计算Hd_奇数
else
Hd(i+1)=0;
end
end
Hd(1,nc+2)=Hd_value*double(limit(Hd_Odd,k,nc+1)); %改善Hd—1
Hd(1,N-nc)=Hd_value*double(limit(Hd_Odd,k,N-nc-1));
% Hd(1,nc+3)=0.01*double(limit(Hd_Odd,k,nc+2)); %改善Hd—2
% Hd(1,N-nc-1)=0.01*double(limit(Hd_Odd,k,N-nc-2));
else %N为偶数
for i=0:N-1
if i<=(N/2-1)
if i<=nc
Hd(i+1)=double(limit(Hd_Even_Front,k,i)); %计算Hd_偶数
else
Hd(i+1)=0;
end
elseif i==N/2
Hd(i+1)=0; %计算Hd_偶数
else
if i>=N-nc
Hd(i+1)=double(limit(Hd_Even_Back,k,i)); %计算Hd_偶数
else
Hd(i+1)=0;
end
end
end
Hd(1,nc+2)=Hd_value*double(limit(Hd_Even_Front,k,nc+1)); %改善Hd—1
Hd(1,N-nc)=Hd_value*double(limit(Hd_Even_Back,k,N-nc-1));
% Hd(1,nc+3)=0.01*double(limit(Hd_Even_Front,k,nc+2)); %改善Hd—2
% Hd(1,N-nc-1)=0.01*double(limit(Hd_Even_Back,k,N-nc-2));
end
hn=ifft(Hd,N); %计算hn
%-----------------------------------------------------------------
%求频率响应
%-------------------
%b=hn; %FIR滤波器的系数
Hw_0=sum(hn); %w=0处的IDTFT(hn)
b=hn/Hw_0; %FIR滤波器的归一化系数
a=1;
[H,w]=freqz(b,a); %返回w为数字频率值
%-----------------------------------------------------------------
%幅频响应&&相频响应
Hr=abs(H); %幅频响应
dBHr=20*log10(Hr); %幅频响应(dB)
%dBHr=Hr;
Hphase=angle(H); %相频响应
Hphase=unwrap(Hphase); % 解卷绕
%-----------------------------------------------------------------
%图形-数据可视化
set(figure(1),'Name','时域实数序列图')
stem(0:M,hn,'.');grid on;
axis([-0.1 M+0.1 -1.1*max(abs(hn)) 1.1*max(abs(hn))]); %限定横坐标和纵坐标的显示范围
title('FIR-频率抽样法-DLPF-时域实数序列');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude value ');
xlabel(' Discrete Time (s)');
%-----------------------------------------------------------------
set(figure(2),'Name','频域序列序列图')
subplot(311)
stem(0:M,real(Hd),'.');grid on; %画频率响应的实部序列--关于中点偶对称
axis([-0.1 M+0.1 -1.1*max(abs(real(Hd))) 1.1*max(abs(real(Hd)))]); %限定横坐标和纵坐标的显示范围
title('FIR-频率抽样法-DLPF-频域值的实部序列（关于中点偶对称）');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude value ');
xlabel(' Discrete Time (s)');
subplot(312)
stem(0:M,imag(Hd),'.');grid on; %画频率响应的虚部序列--关于中点奇对称
axis([-0.1 M+0.1 -1.1*max(abs(imag(Hd))) 1.1*max(abs(imag(Hd)))]); %限定横坐标和纵坐标的显示范围
title('FIR-频率抽样法-DLPF-频域值的虚部序列（关于中点奇对称）');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude value ');
xlabel(' Discrete Time (s)');
subplot(313)
stem(0:M,abs(Hd),'.');grid on; %画频率响应的幅值序列--关于中点偶对称
axis([-0.1 M+0.1 -1.1*max(abs(Hd)) 1.1*max(abs(Hd))]); %限定横坐标和纵坐标的显示范围
title('FIR-频率抽样法-DLPF-频域摸值序列');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude value ');
xlabel(' Discrete Time (s)');
%-----------------------------------------------------------------
set(figure(3),'Name','相频响应图')
plot(w,Hphase,'');grid on; %画相频响应虚线
title('FIR-频率抽样法-DLPF-相频响应');
ylabel(' Phase Freq. Res.');
xlabel(' Digital Freq (rad)');
%-----------------------------------------------------------------
set(figure(4),'Name','幅频响应(dB)图')
%subplot(211)
plot(w,dBHr,'');grid on; %画归一化的幅频率响应曲线
axis([0 pi 1.1*min(dBHr) 1.01*max(dBHr)]); %限定横坐标和纵坐标的显示范围
%axis([0 pi/4 -5 1.01*max(dBHr)]);
title('FIR-频率抽样法-DLPF-幅频响应(dB)');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude Freq. Res.(dB)');
xlabel(' Digital Freq (rad)');