课程编号实验项目序号本科学生实验卡和实验报告信息科学与工程学院通信工程专业2013级1301班课程名称:数字信号处理实验项目:FIR滤波器设计2015~~2016学年第二学期学号:0104_ 姓名:__ _王少丹_ __ 专业年级班级: ____通信1301_________四合院___ 实验室组别________ 实验日期 __2016 年_ 6__ 月___12_ 日(2) 求滤波器的单位抽样响应、频率响应,并绘制波形。
e. 设计FIR 数字带通滤波器,技术指标为:下阻带边缘:ωst1=0.2π,δs1=20dB,下通带边缘:ωp1=0.4π,δp1=1dB;上通带边缘:ωp2=0.6π,δp1=1dB,上阻带边缘:ωst2=0.8π,δs2=20dB;(1) 通过技术指标,选择一种窗函数进行设计;(2) 求滤波器的单位抽样响应、频率响应,并绘制波形。
f. 设计FIR 数字高通滤波器,技术指标为:通带截止频率为ωp=15π/27,阻带截止频率为ωst=11π/27,通带最大衰减为δ1=2.5dB,阻带最小衰减为δ2=55dB。
(1) 通过技术指标,选择一种窗函数进行设计;(2) 求滤波器的单位抽样响应、频率响应,并绘制波形。
系统一:g. 设计FIR 数字高通滤波器,技术指标为:通带截止频率为ωp=0.6π,阻带截止频率为ωst=0.4π,通带最大衰减为δ1=0.25dB,阻带最小衰减为δ2=40dB。
(1) 通过技术指标,选择一种窗函数进行设计;(2) 求滤波器的单位抽样响应、频率响应,并绘制波形。
h. 滤波器的技术指标为:通带截止频率为ωp=0.6π,阻带截止频率为ωst=0.4π,通带最大衰减为δ1=0.25dB,阻带最小衰减为δ2=40dB。
(1) 通过技术指标,选择一种窗函数设计一个具有π/2 相移的FIR 高通滤波器;(2) 求滤波器的单位抽样响应、频率响应,并绘制波形。
i. 设计FIR 数字带阻滤波器,其技术指标为:低端阻带边缘:ωst1=0.4π,δs1=40dB,低端通带边缘:ωp1=0.2π,δp1=1dB;高端通带边缘:ωp2=0.8π,δp1=1dB,高端阻带边缘:ωst2=0.6π,δs2=40dB;(1) 通过技术指标,选择一种窗函数进行设计;(2) 求滤波器的单位抽样响应、频率响应,并绘制波形。
FIR 滤波器的单位抽样响应为h(n)=1/9{},编制MATLAB 程序求系统的频率采样型结构的系数,并画出频率抽样型结构。
m. 一个理想差分器的频率响应为:用长度为21 的汉宁窗设计一个数字FIR 差分器,并绘制其时域和频率的响应波形。
n. 利用汉宁窗设计一个长度为25 的数字希尔伯特变换器,并绘制它的时域和频域的响应波形。
p. FIR 数字低通滤波器的技术指标为:ωp=0.2π,ωst=0.3π,δ1=0.25dB,δ2=50dB。
利用频率采样方法设计FIR 数字滤波器,并绘制滤波器的单位冲激响应、幅度频率响应的波形。
q. 用窗函数法设计一个线性相位的FIR 数字低通滤波器,其技术指标为:ωwp = 0.2*pi;ws = 0.4*pi;tr_width = ws-wp;M = ceil(6.6*pi/tr_width)+1;n = [0:1:M-1];wc = (ws+wp)/2;hd = ideal_lp(wc,M);w_ham = (hamming(M))';h = hd.*w_ham;[db,mag,pha,H,w] = freqz_m3(h,[1]);delta_w = 2*pi/1000;Rp = -(min(db(1:1:wp/delta_w+1)));As = -round(max(db(ws/delta_w+1:1:501)));figure(1)subplot(221)stem(n,hd);title('Ideal Impulse Rresponse')axis([0 M-1 -0.1 0.3]);xlabel('n');ylabel('hd(n)')subplot(222)stem(n,w_ham);title('Hamming Window')axis([0 M-1 0 1.1]);xlabel('n');ylabel('w(n)')subplot(223)stem(n,h);title('Actual Impuse Response')axis([0 M-1 -0.1 0.3]);xlabel('n');ylabel('h(n)')subplot(224)plot(w/pi,db);title('Magnitude Response in db');gridaxis([0 1 -100 10]);xlabel('frequence in \pi unit');ylabel('decibels') %solution of problem (b)Rn = stepseq(0,0,9)n = 0:1:9dw = w(2)-(1)w = [-fliplr(w),w(2:501)]H = [fliplr(H),H(2:501)]xw = Rn*(exp(-j)).^(n'*w)answ = H.*xwRnn = answ*exp(-j).^(w'*n)*dw/(2*pi)figure(2)subplot(221)stem(n,Rn)xlabel('n');ylabel('Rn');title('input rectangle signal in T-domain')subplot(222)stem(w/pi,xw)xlabel('w/\pi');ylabel('Rw');title('input rectangle signal in F-domain') subplot(223)plot(w/pi,abs(answ))xlabel('w/\pi');ylabel('Magnitude of H(w)');title('output signal in F-domain after filteration')subplot(224)stem(n,abs(Rnn))xlabel('n');ylabel('Rnn');title('output signal in T-domain after filteration') M = 25; alpha = (M-1)/2; n = 0:M-1;hd = (2/pi)*((sin((pi/2)*(n-alpha)).^2)./(n-alpha)); hd(alpha+1)=0;w_han = (hann(M))'; h = hd .* w_han; [Hr,w,P,L] = Hr_Type3(h);subplot(2,2,1); stem(n,hd); title('Ideal Impulse Response')axis([-1 M -1.2 1.2]); xlabel('n'); ylabel('hd(n)');subplot(2,2,2); stem(n,w_han);title('Hann Window')axis([-1 M 0 1.2]); xlabel('n'); ylabel('w(n)')subplot(2,2,3); stem(n,h);title('Actual Impulse Response')axis([-1 M -1.2 1.2]); xlabel('n'); ylabel('h(n)')w = w'; Hr = Hr';w = [-fliplr(w), w(2:501)]; Hr = [-fliplr(Hr), Hr(2:501)];subplot(2,2,4);plot(w/pi,Hr); title('Amplitude Response');grid;xlabel('frequency in pi units'); ylabel('Hr'); axis([-1 1 -1.1 1.1]);M = 20; alpha = (M-1)/2; l = 0:M-1; wl = (2*pi/M)*l;Hrs = [1,1,1,zeros(1,15),1,1]; %Ideal Amp Res sampledHdr = [1,1,0,0]; wdl = [0,0.25,0.25,1]; %Ideal Amp Res for plottingk1 = 0:floor((M-1)/2); k2 = floor((M-1)/2)+1:M-1;angH = [-alpha*(2*pi)/M*k1, alpha*(2*pi)/M*(M-k2)];H = Hrs.*exp(j*angH); h = real(ifft(H,M));[db,mag,pha,grd,w] = freqz_m(h,1); [Hr,ww,a,L] = Hr_Type2(h);subplot(2,2,1);plot(wl(1:11)/pi,Hrs(1:11),'o',wdl,Hdr);axis([0,1,-0.1,1.1]); title('Frequency Samples: M=20')xlabel('frequency in pi units'); ylabel('Hr(k)');subplot(2,2,2); stem(l,h); axis([-1,M,-0.1,0.3])title('Impulse Response'); xlabel('n'); ylabel('h(n)');subplot(2,2,3); plot(ww/pi,Hr,wl(1:11)/pi,Hrs(1:11),'o');axis([0,1,-0.2,1.2]); title('Amplitude Response')xlabel('frequency in pi units'); ylabel('Hr(w)')subplot(2,2,4);plot(w/pi,db); axis([0,1,-60,10]); gridtitle('Magnitude Response'); xlabel('frequency in pi units');ylabel('Decibels');Wpl=0.2*pi;Wph=0.8*pi;Wsl=0.4*pi;Wsh=0.6*pi;tr_width=min((Wsl-Wpl),(Wph-Wsh));M=ceil(6.2*pi/tr_width)n=0:1:M-1;Wcl=(Wsl+Wpl)/2;Wch=(Wsh+Wph)/2;hd=ideal_bs(Wcl,Wch,M);w_ham=(hanning(M))';h=hd.*w_ham;[db,mag,pha,w]=freqz_m2(h,[1]);delta_w=2*pi/1000;Ap=-(min(db(1:1:Wpl/delta_w+1)))As=-round(max(db(Wsl/delta_w+1:1:Wsh/delta_w+1)))subplot(221)stem(n,hd);title('Ideal Impulse Rresponse')axis([0 M-1 -0.1 0.7]);xlabel('n');ylabel('hd(n)')subplot(222)stem(n,w_ham);title('Hamming Window')axis([0 M-1 0 1.1]);xlabel('n');ylabel('w(n)')subplot(223)stem(n,h);title('Actual Impuse Response')axis([0 M-1 -0.1 0.7]);xlabel('n');ylabel('h(n)')subplot(224)plot(w/pi,db);title('Magnitude Response in db');gridaxis([0 1 -100 10]);xlabel('frequence in \pi unit');ylabel('decibels'); clear allfigureWpl=0.2*pi;Wph=0.8*pi;Wsl=0.4*pi;Wsh=0.6*pi;tr_width=min((Wsl-Wpl),(Wph-Wsh));M=ceil(6.2*pi/tr_width)n=0:1:M-1;Wcl=(Wsl+Wpl)/2;Wch=(Wsh+Wph)/2;hd=ideal_bs(Wcl,Wch,M);w_ham=(hanning(M))';h=hd.*w_ham;[db,mag,pha,w]=freqz_m2(h,[1]);delta_w=2*pi/1000;Ap=-(min(db(1:1:Wpl/delta_w+1)))As=-round(max(db(Wsl/delta_w+1:1:Wsh/delta_w+1)))subplot(221)stem(n,hd);title('Ideal Impulse Rresponse')axis([0 M-1 -0.1 0.7]);xlabel('n');ylabel('hd(n)')subplot(222)stem(n,w_ham);title('Hamming Window')axis([0 M-1 0 1.1]);xlabel('n');ylabel('w(n)')subplot(223)stem(n,h);title('Actual Impuse Response')axis([0 M-1 -0.1 0.7]);xlabel('n');ylabel('h(n)')subplot(224)plot(w/pi,db);title('Magnitude Response in db');gridaxis([0 1 -100 10]);xlabel('frequence in \pi unit');ylabel('decibels'); figurews1=0.2*pi;wp1=0.35*pi;ws2=0.8*pi;wp2=0.65*pi;Ap=60;Rp=1;tr_width=min((wp1-ws1),(ws2-wp2));M=ceil(11*pi/tr_width);n=[0:1:M-1];wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);w_bla=(blackman(M))';h=hd.*w_bla;[H,W]=freqz(h,1);subplot(2,2,1);stem(n,hd);title('ÀíÏëÂö³å³éÑù');subplot(2,2,2);stem(n,w_bla);title('²¼À³¿ËÂü´°');subplot(2,2,3);stem(n,h);title('ʵ¼ÊÂö³å³éÑù');subplot(2,2,4);plot(W/pi,20*log10(abs(H)));title('·ù¶ÈÏìÓ¦£¨db£©'); clear allfigurews1=0.2*pi;wp1=0.4*pi;ws2=0.8*pi;wp2=0.6*pi;Ap=60;Rp=1;tr_width=min((wp1-ws1),(ws2-wp2));M=ceil(11*pi/tr_width);n=[0:1:M-1];wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);w_bla=( hamming(M))';h=hd.*w_bla;[H,W]=freqz(h,1);subplot(2,2,1);stem(n,hd);title('ÀíÏëÂö³å³éÑù');subplot(2,2,2);stem(n,w_bla);title('º£Ã÷´°');subplot(2,2,3);stem(n,h);title('ʵ¼ÊÂö³å³éÑù');subplot(2,2,4);plot(W/pi,20*log10(abs(H)));title('·ù¶ÈÏìÓ¦£¨db£©'); clear allfigurews1=0.2*pi;wp1=0.4*pi;ws2=0.8*pi;wp2=0.6*pi;Ap=20;Rp=1;tr_width=min((wp1-ws1),(ws2-wp2));M=ceil(11*pi/tr_width);n=[0:1:M-1];wc1=(ws1+wp1)/2;wc2=(wp2+ws2)/2;hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);w_bla=( boxcar(M))';h=hd.*w_bla;[H,W]=freqz(h,1);subplot(2,2,1);stem(n,hd);title('ÀíÏëÂö³å³éÑù');subplot(2,2,2);stem(n,w_bla);title('¾ØÐδ°');subplot(2,2,3);stem(n,h);title('ʵ¼ÊÂö³å³éÑù');subplot(2,2,4);plot(W/pi,20*log10(abs(H)));title('·ù¶ÈÏìÓ¦£¨db£©'); clear allfigureAs=55;ws=11*pi/27;wp=15*pi/27;tr_width=wp-ws;%¼ÆËã¹ý¶É´øM=ceil((As-7.95)*2*pi/(14.36*tr_width)+1)+1; %°´¿-Ôó´°¼ÆËãÂ˲¨Æ÷µÄ³¤¶Èdisp(['Â˲¨Æ÷µÄ³¤¶È',num2str(M)]);beta=0.1102*(As-8.7); %¼ÆËã¿-Ôó´°µÄbetaÖµn=[0:1:M-1];disp(['ÏßÐÔÏàλÂ˲¨Æ÷',num2str(beta)]);w_kai=(kaiser(M,beta))';%Çó¿-Ôó´°º¯Êýwc=(ws+wp)/2;hd=ideal_lp(pi,M)-ideal_lp(wc,M); %ÇóÀíÏëÂö³åÏìÓ¦h=hd.*w_kai;[db,mag,pha,grd,w]=freqz_m(h,[1]);delta_w=2*pi/1000;Rp=-(min(db(wp/delta_w+1:1:501)));disp(['ʵ¼Êͨ´ø²¨¶¯Îª',num2str(Rp)]);As=-round(max(db(1:1:ws/delta_w+1)));disp(['×îС×è´øË¥¼õΪ?',num2str(As)]);subplot(2,2,1);stem(n,hd);title('ÀíÏëÂö³åÏìÓ¦');axis([0 M-1 -0.4 0.8]);ylabel('hd(n)');subplot(2,2,2);stem(n,w_kai);title('¿-Ôó´°');axis([0 M-1 0 1.1]);ylabel('wd(n)');subplot(2,2,3);stem(n,h);title('ʵ¼ÊÂö³åÏìÓ¦');axis([0 M-1 -0.4 0.8]);xlabel('n');ylabel('h(n)');subplot(2,2,4);plot(w/pi,db);title('?·ù¶ÈÏìÓ¦/dB');axis([0 1 -100 10]);grid;xlabel('ÒÔpiΪµ¥Î»µÄƵÂÊ');ylabel('?·Ö±´Êý/dB');clear allfigure%ÔÚmatlabÖÐн¨º¯Êýg£¬´úÂëÈçÏ£ºWp=0.6*pi;Ws=0.4*pi;tr_width=Wp-Ws;M=ceil(6.2*pi/tr_width);n=0:1:M-1;Wc=(Ws+Wp)/2;hd=ideal_lp(pi,M)-ideal_lp(Wc,M);w_ham=(hanning(M))';h=hd.*w_ham;[db,mag,pha,w]=freqz_m2(h,[1]);delta_w=2*pi/1000;Ap=-(min(db(Wp/delta_w+1:1:501)))As=-round(max(db(1:1:Ws/delta_w+1)))subplot(221)stem(n,hd);title('Ideal Impulse Rresponse')axis([0 M-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)') subplot(222)stem(n,w_ham);title('Hamming Window')axis([0 M-1 0 1.1]);xlabel('n');ylabel('w(n)') subplot(223)stem(n,h);title('Actual Impuse Response')axis([0 M-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')subplot(224)plot(w/pi,db);title('Magnitude Response in db');gridaxis([0 1 -100 10]);xlabel('frequence in \pi unit');ylabel('decibels'); clear allfigureWp=0.6*pi;Ws=0.4*pi;tr_width=Wp-Ws;M=ceil(6.2*pi/tr_width)n=0:1:M-1;Wc=(Ws+Wp)/2;hd=ideal_lp(pi,M)-ideal_lp(Wc,M);w_ham=(hanning(M));h=hd.*rot90(w_ham);[db,mag,pha,w]=freqz_m2(h,[1]);delta_w=2*pi/1000;Ap=-(min(db(Wp/delta_w+1:1:501)))As=-round(max(db(1:1:Ws/delta_w+1)))figuresubplot(221)stem(n,hd);title('Ideal Impulse Rresponse')axis([0 M-1 -0.4 0.5]);xlabel('n');ylabel('hd(n)')subplot(222)stem(n,w_ham);title('Hamming Window')axis([0 M-1 0 1.1]);xlabel('n');ylabel('w(n)')subplot(223)stem(n,h);title('Actual Impuse Response')axis([0 M-1 -0.4 0.5]);xlabel('n');ylabel('h(n)')subplot(224)plot(w/pi,db);title('Magnitude Response in db');gridaxis([0 1 -100 10]);xlabel('frequence in \pi unit');ylabel('decibels'); %ÓÃHamming´°º¯ÊýÉè¼ÆFIRÊý×ÖÂ˲¨Æ÷figurewp=0.2*pi;ws=0.3*piN=61n=[0:1:N-1]wc=(ws+wp)/2;%ÀíÏëµÍͨÂ˲¨Æ÷hd=ideal_lp(wc,N);%ÀíÏëµÍͨµÄ³å¼¤ÏìÓ¦w_ham=(hamming(N))';h=hd.*w_ham;%FIRÂ˲¨Æ÷³å¼¤ÏìÓ¦[db,mag,pha,grd,w]=freqz_m(h,[1]);delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1)))%ʵ¼ÊµÄͨ´øË¥¼õAs=-round(max(db(ws/delta_w+1:1:501)))%ʵ¼ÊµÄ×îС×è´øË¥¼õsubplot(221);stem(n,hd);title('ÀíÏë³å¼¤ÏìÓ¦')axis([0 N-1 -0.1 0.3]);xlabel('n');ylabel('hd(n)')subplot(222);stem(n,w_ham);title('hamming´°')axis([0 N-1 0 1.1]);xlabel('n');ylabel('w(n)');subplot(223);stem(n,h);title('ʵ¼Ê³å¼¤ÏìÓ¦')axis([0 N-1 -0.1 0.3]);xlabel('n');ylabel('h(n)')subplot(224);plot(w/pi,db);axis([0 0.8 -100 0]);xlabel('ÒÔPIΪµ¥Î»µÄƵÂÊ');ylabel('¶ÔÊý·ù¶È/db'); clear allfigure%ÓÃKaiser´°º¯ÊýÉè¼ÆFIRÊý×ÖÂ˲¨Æ÷wp=0.2*pi;ws=0.3*pi;As=50tr_width=ws-wpN=ceil((As-7.95)/(14.36*tr_width/(2*pi))+1)+1n=[0:1:N-1]beta=0.1102*(As-8.7)wc=(wp+ws)/2%ÀíÏëµÍͨµÄ½ØÖ¹ÆµÂÊhd=ideal_lp(wc,N)w_kai=(kaiser(N,beta))'h=hd.*w_kai[db,mag,pha,grd,w]=freqz_m(h,[1])delta_w=2*pi/1000Rp=-(min(db(1:1:wp/delta_w+1)))%?ʵ¼ÊµÄͨ´øË¥¼õAs=-round(max(db(ws/delta_w+1:1:501)))%ʵ¼ÊµÄ×îС×è´øË¥¼õsubplot(211);plot(w/pi,db);title('¿-É-´°·ù¶ÈÏìÓ¦(dB)');gridaxis([0 0.5 -100 0])ylabel('¶ÔÊý·ù¶È/db');xlabel('ÒÔ/piΪµ¥Î»µÄƵÂÊ')subplot(212);plot(w/pi,pha);title('ÏàλÏìÓ¦');gridaxis([0 0.5 -4 4])ylabel('Ïàλ');xlabel('ÒÔ/piΪµ¥Î»µÄƵÂÊ')figureh1=[-4,1,-1,-2,5,6,5,-2,-1,1,-4];M=length(h1);n=0:M-1;[A,w,type,tao]=amplres(h1);type subplot(2,1,1),stem(n,h1);title('³å¼¤ÏìÓ¦h1');ylabel('h(n)');xlabel('n');subplot(2,1,2),plot(w/pi,A);ylabel('A');xlabel('\pi');title('?·ùƵÏìÓ¦');figurerz=roots(h1)for i=1:8r(i)=1/rz(i);endr'zplane(h1,1);title('h1Á㼫µãͼ');clear allfigureh2=[-4,1,-1,-2,5,6,6,5,-2,-1,1,-4]; M=length(h2);n=0:M-1;[A,w,type,tao]=amplres(h2);type subplot(2,1,1),stem(n,h2);title('³å¼¤ÏìÓ¦h2');ylabel('h(n)');xlabel('n');subplot(2,1,2),plot(w/pi,A);ylabel('A');xlabel('\pi');title('·ùƵÏìÓ¦');figurerz=roots(h2)for i=1:8r(i)=1/rz(i);endr'zplane(h2,1);title('?h2Á㼫µãͼ');clear allfigureh3=[-4,1,-1,-2,5,0,-5,2,1,-1,4];M=length(h3);n=0:M-1;[A,w,type,tao]=amplres(h3);type subplot(2,1,1),stem(n,h3);title('³å¼¤ÏìÓ¦h3');ylabel('h(n)');xlabel('n');subplot(2,1,2),plot(w/pi,A);ylabel('A');xlabel('\pi');title('·ùƵÏìÓ¦?');figurerz=roots(h3)for i=1:8r(i)=1/rz(i);endr'zplane(h3,1);title('h3Á㼫µãͼ');clear allfigureh4=[-4,1,-1,-2,5,6,-6,-5,2,1,-1,4]; M=length(h4);n=0:M-1;[A,w,type,tao]=amplres(h4);type subplot(2,1,1),stem(n,h4);title('³å¼¤ÏìÓ¦h4');ylabel('h(n)');xlabel('n');subplot(2,1,2),plot(w/pi,A);ylabel('A');xlabel('\pi');title('?·ùƵÏìÓ¦?');figurerz=roots(h4)for i=1:8r(i)=1/rz(i);endr'zplane(h4,1);title('h4Á㼫µãͼ');h = [1,2,3,2,1];M = 5;n = (M-3)/2;a = 2*h(3-n)k = -500:500;K = 500;w = k*pi/K;p = a*cos(w*n);plot(w,p);ylabel('p(w)'); xlabel('w');function[db,mag,pha,H,w]=freqz_m3(b,a)[H,w]=freqz(b,a,1000,'whole');H=(H(1:1:501))'; w=(w(1:1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);endfunction [db,mag,pha,grd,w] = freqz_m(b,a);[H,w] = freqz(b,a,1000,'whole');H = (H(1:1:501))'; w = (w(1:1:501))';mag = abs(H);db = 20*log10((mag+eps)/max(mag));pha = angle(H);grd = grpdelay(b,a,w);function [db, mag, pha, w]=freqz_m2(b,a);%Modified version of freqz subroutine[H,w]=freqz(b,a,1000,'whole');H=(H(1:501))'; w=(w(1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);function [Hr,w,b,L] = Hr_Type2(h);M = length(h); L = M/2;b = 2*[h(L:-1:1)]; n = [1:1:L]; n = n-0.5;w = [0:1:500]'*pi/500; Hr = cos(w*n)*b';endfunction [Hr,w,c,L] = Hr_Type3(h);% Computes Amplitude response Hr(w) of a Type-3 LP FIR filter % -----------------------------------------------------------% [Hr,w,c,L] = Hr_Type3(h)% Hr = Amplitude Response% w = frequencies between [0 pi] over which Hr is computed % c = Type-3 LP filter coefficients% L = Order of Hr% h = Type-3 LP impulse response%M = length(h); L = (M-1)/2;c = [2*h(L+1:-1:1)]; n = [0:1:L];w = [0:1:500]'*pi/500; Hr = sin(w*n)*c';end。