数字信号处理研究学习.doc

数字信号处理课程研究性学习报告DFT近似计算信号频谱专题研讨【目的】(1) 掌握利用DFT近似计算不同类型信号频谱的原理和方法。(2) 理解误差产生的原因及减小误差的方法。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目1】 基本题1利用DFT分析x(t)=Acos(2f1t)+Bcos(2f2t)的频谱，其中f1=100Hz，f2=120Hz(1)A=B=1; (2)A=1,B=0.2。要求选择不同的窗函数。【题目分析】由取样定理知，要使信号频谱不混叠，则抽样频率不小于最高频率的两倍，即。而要满足信号分辨率的要求，抽样点数。FFT点数要求.以下用控制变量的办法，分别研究fs、N和L对频谱误差的影响。（1）研究抽样频率fs对频谱误差的影响时，固定N=100(,此处取N=100可认为N带来的误差可忽略)。【仿真程序1】N=100;L=512;35T=1/50;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(311);plot(f,abs(y1),linewidth,1.5);title(图表一抽样频率fs取值原则研究, fsfmax );xlabel(f);ylabel(幅度值);N=100;L=512;T=1/125;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(312);plot(f,abs(y1),linewidth,1.5);title(图表一抽样频率fs取值原则研究, fmaxfs2*fmax );xlabel(f);ylabel(幅度值);由结果可知：在fsfmax 频率混叠 ，fmaxfs2*fmax时均无混叠，这是因为余弦信号是带限信号，只要抽样频率大于最高频率便无混叠。（2）研究时域抽样点N对频谱误差的影响时，固定fs=600（可保证无混叠）。改变N：【仿真程序2】N=10;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);figure(2);subplot(411);plot(f,abs(y1),linewidth,1.5);title(图表二抽样点数N取值原则研究, Nfs/delta f,N=10);xlabel(f);ylabel(幅度值);N=20;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(412);plot(f,abs(y1),linewidth,1.5);title(图表二抽样点数N取值原则研究, Nfs/delta f,N=20);xlabel(f);ylabel(幅度值);N=30;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(413);plot(f,abs(y1),linewidth,1.5);title(图表二抽样点数N取值原则研究, N=fs/delta f,N=30);xlabel(f);ylabel(幅度值);N=60;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(414);plot(f,abs(y1),linewidth,1.5);title(图表二抽样点数N取值原则研究, N=2*fs/delta f,N=60);xlabel(f);ylabel(幅度值);可见，当时，虽然两谱峰有混叠，但基本上可将其分辨出，而当时，两者完全分离，无混叠。（3）研究FFT点数L对频谱误差的影响时，固定N=60,fs=600；改变L:【仿真程序3】N=60;L=30;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);figure(3);subplot(411);plot(f,abs(y1),linewidth,1.5);title(图表3DFT点数L取值原则研究, LN L=62);xlabel(f);ylabel(幅度值);N=60;L=128;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(413);plot(f,abs(y1),linewidth,1.5);title(图表3DFT点数L取值原则研究, LN L=128);xlabel(f);ylabel(幅度值);N=60;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+cos(2*pi*f2*t);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(414);plot(f,abs(y1),linewidth,1.5);title(图表3DFT点数L取值原则研究, LN L=512);xlabel(f);ylabel(幅度值);由实验结果可知，即便L=N,或者稍大于N,也不能保证所得到的曲线是光滑的，在保证运算速度的基础上，应尽量增大FFT点数以确保频谱信息完整。2.而当A=1;B=0.2时由于混叠误差、栅栏效应与两者等值时是一样的道理，因而此处只研究窗函数的长度带来的误差。改变N值：10 20 30 60 90； 【仿真程序4】N=10;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);x=x.*hamming(N);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);figure(2);subplot(511);plot(f,abs(y1),linewidth,1.5);title(图表二抽样点数N取值原则研究, Nfs/delta f,N=10);xlabel(f);ylabel(幅度值);N=20;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);x=x.*hamming(N);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(512);plot(f,abs(y1),linewidth,1.5);title(图表二抽样点数N取值原则研究, Nfs/delta f,N=20);xlabel(f);ylabel(幅度值);N=30;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);x=x.*hamming(N);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(513);plot(f,abs(y1),linewidth,1.5);title(图表二抽样点数N取值原则研究, N=fs/delta f,N=30);xlabel(f);ylabel(幅度值);N=60;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);x=x.*hamming(N);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(514);plot(f,abs(y1),linewidth,1.5);title(图表二抽样点数N取值原则研究, N=2*fs/delta f,N=60);xlabel(f);ylabel(幅度值);N=90;L=512;T=1/600;t=(0:N-1)*T;f1=100;f2=120;x=cos(2*pi*f1*t)+0.2*cos(2*pi*f2*t);x=x.*hamming(N);y1=fftshift(fft(x,L);ws=2*pi/T;f=(-0.5*ws+(0:L-1)*ws/L)/(2*pi);subplot(515);plot(f,abs(y1),linewidth,1.5);title(图表二抽样点数N取值原则研究, N=2*fs/delta f,N=90);xlabel(f);ylabel(幅度值);可见使用汉明窗时也未必就能分得开，这只是一个理论上的分析值，究竟能否分开还需进行实验。实验结论：综上，连续时间信号选择合适DFT参数原则为为：（1）（非带限信号）；（2）窗函数的长度参考值是：矩形窗，汉明窗，（但究竟取多少还需进行实验而确定）； （3）（在保证运算速度的基础上，应尽量增大L以确保频谱信息完整。）【研讨题目2】 基本题2已知一离散序列为 (1)用L=32点DFT计算该序列的频谱，求出频谱中谱峰的频率；(2)对序列进行补零，然后分别用L=64、128、256、512点DFT计算该序列的频谱，求出频谱中谱峰的频率；(3)讨论所获得的结果，给出你的结论。该结论对序列的频谱计算有何指导意义？【题目分析】(1)(2)两问要求使用不同的DFT点数进行频谱分析，并且分析谱峰值与理论值的差别，进而分析DFT点数对频谱计算的影响。其程序与结果如下.【仿真程序】N=32;L=32;k=0:N-1;x=sin(0.2*pi*k);y=fft(x,L);w=2*(0:L-1)/L;subplot(511);plot(w,abs(y),linewidth,2);title(L=32);xlabel(omega/pi);ylabel(频谱幅度值);L=64;k=0:N-1;x=sin(0.2*pi*k);y=fft(x,L);w=2*(0:L-1)/L;subplot(512);plot(w,abs(y),linewidth,2);xlabel(omega/pi);ylabel(频谱幅度值);title(L=64);L=128;k=0:N-1;x=sin(0.2*pi*k);y=fft(x,L);w=2*(0:L-1)/L;subplot(513);plot(w,abs(y),linewidth,2);xlabel(omega/pi);ylabel(频谱幅度值);title(L=128);L=256;k=0:N-1;x=sin(0.2*pi*k);y=fft(x,L);w=2*(0:L-1)/L;subplot(514);plot(w,abs(y),linewidth,2);xlabel(omega/pi);ylabel(频谱幅度值);title(L=256);L=512;k=0:N-1;x=sin(0.2*pi*k);y=fft(x,L);w=2*(0:L-1)/L;subplot(515);plot(w,abs(y),linewidth,2);xlabel(omega/pi);ylabel(频谱幅度值);title(L=512);峰值对应频率的计算：由于离散非周期序列的L点DFT是对其在（0，2）上的L点等间隔抽样（1）L=32；由此计算抽样后间隔为，可知抽样值为0.0625*k，即0 0.0625 0. 0.1250 0.1875 0.2500 0.3125 0.3750 0.4375 0.5000 0.5625 0.6250 0.6875 0.7500 0.8125 0.8750 0.9375 1.0000 1.0625 1.1250 1.1875 1.2500 1.3125 1.3750 1.4375 1.5000 1.5625 1.6250 1.6875 1.7500 1.8125 1.8750 1.9375比对仿真所得图表可知(即最接近0.2及1.8)，峰值对应的频率为0.1875，1.8125。（2）L=64；由此计算抽样后间隔为，所以抽样值为0.03125*k，即 0 0.0313 0.0625 0.0938 0.1250 0.1563 0.1875 0.2188 0.25000.2813 0.3125 0.3438 0.3750 0.4063 0.4375 0.4688 0.5000 0.53130.5625 0.5938 0.6250 0.6563 0.6875 0.7188 0.7500 0.7813 0.81250.8438 0.8750 0.9063 0.9375 0.9688 1.0000 1.0313 1.0625 1.09381.1250 1.1563 1.1875 1.2188 1.2500 1.2813 1.3125 1.3438 1.37501.4063 1.4375 1.4688 1.5000 1.5313 1.5625 1.5938 1.6250 1.65631.6875 1.7188 1.7500 1.7813 1.8125 1.8438 1.8750 1.9063 1.93751.9688比对仿真所得图表可知，峰值对应的频率为0.1875，1.8125。（3）L=128；由此计算抽样后间隔为，所以抽样值为1/64*k，即 0 0.0156 0.0313 0.0469 0.0625 0.0781 0.0938 0.1094 0.1250 0.1406 0.1563 0.1719 0.1875 0.2031 0.2188 0.2344 0.2500 0.2656 0.2813 0.2969 0.3125 0.3281 0.3438 0.3594 0.3750 0.3906 0.4063 0.4219 0.4375 0.4531 0.4688 0.4844 0.5000 0.5156 0.5313 0.5469 0.5625 0.5781 0.5938 0.6094 0.6250 0.6406 0.6563 0.6719 0.6875 0.7031 0.7188 0.7344 0.7500 0.7656 0.7813 0.7969 0.8125 0.8281 0.8438 0.8594 0.8750 0.8906 0.9063 0.9219 0.9375 0.9531 0.9688 0.9844 1.0000 1.0156 1.0313 1.0469 1.0625 1.0781 1.0938 1.1094 1.1250 1.1406 1.1563 1.1719 1.1875 1.2031 1.2188 1.2344 1.2500 1.2656 1.2813 1.2969 1.3125 1.3281 1.3438 1.3594 1.3750 1.3906 1.4063 1.4219 1.4375 1.4531 1.4688 1.4844 1.5000 1.5156 1.5313 1.5469 1.5625 1.5781 1.5938 1.6094 1.6250 1.6406 1.6563 1.6719 1.6875 1.7031 1.7188 1.7344 1.7500 1.7656 1.7813 1.7969 1.8125 1.8281 1.8438 1.8594 1.8750 1.8906 1.9063 1.9219 1.9375 1.9531 1.9688 1.9844比对仿真所得图表可知，峰值对应的频率为0.2031，1.7969。（4）L=256；由此计算抽样后间隔为，所以抽样值为1/128*k，即 0 0.0078 0.0156 0.0234 0.0313 0.0391 0.0469 0.0547 0.0625 0.0703 0.0781 0.0859 0.0938 0.1016 0.1094 0.1172 0.1250 0.1328 0.1406 0.1484 0.1563 0.1641 0.1719 0.1797 0.1875 0.1953 0.2031 0.2109 0.2188 0.2266 0.2344 0.2422 0.2500 0.2578 0.2656 0.2734 0.2813 0.2891 0.2969 0.3047 0.3125 0.3203 0.3281 0.3359 0.3438 0.3516 0.3594 0.3672 0.3750 0.3828 0.3906 0.3984 0.4063 0.4141 0.4219 0.4297 0.4375 0.4453 0.4531 0.4609 0.4688 0.4766 0.4844 0.4922 0.5000 0.5078 0.5156 0.5234 0.5313 0.5391 0.5469 0.5547 0.5625 0.5703 0.5781 0.5859 0.5938 0.6016 0.6094 0.6172 0.6250 0.6328 0.6406 0.6484 0.6563 0.6641 0.6719 0.6797 0.6875 0.6953 0.7031 0.7109 0.7188 0.7266 0.7344 0.7422 0.7500 0.7578 0.7656 0.7734 0.7813 0.7891 0.7969 0.8047 0.8125 0.8203 0.8281 0.8359 0.8438 0.8516 0.8594 0.8672 0.8750 0.8828 0.8906 0.8984 0.9063 0.9141 0.9219 0.9297 0.9375 0.9453 0.9531 0.9609 0.9688 0.9766 0.9844 0.9922 1.0000 1.0078 1.0156 1.0234 1.0313 1.0391 1.0469 1.0547 1.0625 1.0703 1.0781 1.0859 1.0938 1.1016 1.1094 1.1172 1.1250 1.1328 1.1406 1.1484 1.1563 1.1641 1.1719 1.1797 1.1875 1.1953 1.2031 1.2109 1.2188 1.2266 1.2344 1.2422 1.2500 1.2578 1.2656 1.2734 1.2813 1.2891 1.2969 1.3047 1.3125 1.3203 1.3281 1.3359 1.3438 1.3516 1.3594 1.3672 1.3750 1.3828 1.3906 1.3984 1.4063 1.4141 1.4219 1.4297 1.4375 1.4453 1.4531 1.4609 1.4688 1.4766 1.4844 1.4922 1.5000 1.5078 1.5156 1.5234 1.5313 1.5391 1.5469 1.5547 1.5625 1.5703 1.5781 1.5859 1.5938 1.6016 1.6094 1.6172 1.6250 1.6328 1.6406 1.6484 1.6563 1.6641 1.6719 1.6797 1.6875 1.6953 1.7031 1.7109 1.7188 1.7266 1.7344 1.7422 1.7500 1.7578 1.7656 1.7734 1.7813 1.7891 1.7969 1.8047 1.8125 1.8203 1.8281 1.8359 1.8438 1.8516 1.8594 1.8672 1.8750 1.8828 1.8906 1.8984 1.9063 1.9141 1.9219 1.9297 1.9375 1.9453 1.9531 1.96091.9688 1.9766 1.9844 1.9922比对仿真所得图表可知，峰值对应的频率为0.2031，1.7969。（5）L=512；由此计算抽样后间隔为，所以抽样值为1/256*k，即 0 0.0039 0.0078 0.0117 0.0156 0.0195 0.0234 0.0273 0.03130.0352 0.0391 0.0430 0.0469 0.0508 0.0547 0.0586 0.0625 0.06640.0703 0.0742 0.0781 0.0820 0.0859 0.0898 0.0938 0.0977 0.10160.1055 0.1094 0.1133 0.1172 0.1211 0.1250 0.1289 0.1328 0.13670.1406 0.1445 0.1484 0.1523 0.1563 0.1602 0.1641 0.1680 0.17190.1758 0.1797 0.1836 0.1875 0.1914 0.1953 0.1992 0.2031 0.20700.2109 0.2148 0.2188 0.2227 0.2266 0.2305 0.2344 0.2383 0.24220.2461 0.2500 0.2539 0.2578 0.2617 0.2656 0.2695 0.2734 0.27730.2813 0.2852 0.2891 0.2930 0.2969 0.3008 0.3047 0.3086 0.31250.3164 0.3203 0.3242 0.3281 0.3320 0.3359 0.3398 0.3438 0.34770.3516 0.3555 0.3594 0.3633 0.3672 0.3711 0.3750 0.3789 0.38280.3867 0.3906 0.3945 0.3984 0.4023 0.4063 0.4102 0.4141 0.41800.4219 0.4258 0.4297 0.4336 0.4375 0.4414 0.4453 0.4492 0.45310.4570 0.4609 0.4648 0.4688 0.4727 0.4766 0.4805 0.4844 0.48830.4922 0.4961 0.5000 0.5039 0.5078 0.5117 0.5156 0.5195 0.52340.5273 0.5313 0.5352 0.5391 0.5430 0.5469 0.5508 0.5547 0.55860.5625 0.5664 0.5703 0.5742 0.5781 0.5820 0.5859 0.5898 0.59380.5977 0.6016 0.6055 0.6094 0.6133 0.6172 0.6211 0.6250 0.62890.6328 0.6367 0.6406 0.6445 0.6484 0.6523 0.6563 0.6602 0.66410.6680 0.6719 0.6758 0.6797 0.6836 0.6875 0.6914 0.6953 0.69920.7031 0.7070 0.7109 0.7148 0.7188 0.7227 0.7266 0.7305 0.73440.7383 0.7422 0.7461 0.7500 0.7539 0.7578 0.7617 0.7656 0.76950.7734 0.7773 0.7813 0.7852 0.7891 0.7930 0.7969 0.8008 0.80470.8086 0.8125 0.8164 0.8203 0.8242 0.8281 0.8320 0.8359 0.83980.8438 0.8477 0.8516 0.8555 0.8594 0.8633 0.8672 0.8711 0.87500.8789 0.8828 0.8867 0.8906 0.8945 0.8984 0.9023 0.9063 0.91020.9141 0.9180 0.9219 0.9258 0.9297 0.9336 0.9375 0.9414 0.94530.9492 0.9531 0.9570 0.9609 0.9648 0.9688 0.9727 0.9766 0.98050.9844 0.9883 0.9922 0.9961 1.0000 1.0039 1.0078 1.0117 1.01561.0195 1.0234 1.0273 1.0313 1.0352 1.0391 1.0430 1.0469 1.05081.0547 1.0586 1.0625 1.0664 1.0703 1.0742 1.0781 1.0820 1.08591.0898 1.0938 1.0977 1.1016 1.1055 1.1094 1.1133 1.1172 1.12111.1250 1.1289 1.1328 1.1367 1.1406 1.1445 1.1484 1.1523 1.15631.1602 1.1641 1.1680 1.1719 1.1758 1.1797 1.1836 1.1875 1.19141.1953 1.1992 1.2031 1.2070 1.2109 1.2148 1.2188 1.2227 1.22661.2305 1.2344 1.2383 1.2422 1.2461 1.2500 1.2539 1.2578 1.26171.2656 1.2695 1.2734 1.2773 1.2813 1.2852 1.2891 1.2930 1.29691.3008 1.3047 1.3086 1.3125 1.3164 1.3203 1.3242 1.3281 1.33201.3359 1.3398 1.3438 1.3477 1.3516 1.3555 1.3594 1.3633 1.36721.3711 1.3750 1.3789 1.3828 1.3867 1.3906 1.3945 1.3984 1.40231.4063 1.4102 1.4141 1.4180 1.4219 1.4258 1.4297 1.4336 1.43751.4414 1.4453 1.4492 1.4531 1.4570 1.4609 1.4648 1.4688 1.47271.4766 1.4805 1.4844 1.4883 1.4922 1.4961 1.5000 1.5039 1.50781.5117 1.5156 1.5195 1.5234 1.5273 1.5313 1.5352 1.5391 1.54301.5469 1.5508 1.5547 1.5586 1.5625 1.5664 1.5703 1.5742 1.57811.5820 1.5859 1.5898 1.5938 1.5977 1.6016 1.6055 1.6094 1.61331.6172 1.6211 1.6250 1.6289 1.6328 1.6367 1.6406 1.6445 1.64841.6523 1.6563 1.6602 1.6641 1.6680 1.6719 1.6758 1.6797 1.68361.6875 1.6914 1.6953 1.6992 1.7031 1.7070 1.7109 1.7148 1.71881.7227 1.7266 1.7305 1.7344 1.7383 1.7422 1.7461 1.7500 1.75391.7578 1.7617 1.7656 1.7695 1.7734 1.7773 1.7813 1.7852 1.78911.7930 1.7969 1.8008 1.8047 1.8086 1.8125 1.8164 1.8203 1.82421.8281 1.8320 1.8359 1.8398 1.8438 1.8477 1.8516 1.8555 1.85941.8633 1.8672 1.8711 1.8750 1.8789 1.8828 1.8867 1.8906 1.89451.8984 1.9023 1.9063 1.9102 1.9141 1.9180 1.9219 1.9258 1.92971.9336 1.9375 1.9414 1.9453 1.9492 1.9531 1.9570 1.9609 1.96481.9688 1.9727 1.9766 1.9805 1.9844 1.9883 1.9922 1.9961比对仿真所得图表可知，峰值对应的频率为0.1992 1.8008。结果分析与结论比较可知，DFT点数越多，谱峰对应的频率值与理论值差别越小，但是当DFT点数较多时，这种差别是很微小的。离散非周期序列的DFT点数越多，计算所得的DFT包含的信息越多，越接近，但是DFT点数过多会增大计算量，减慢运算速度。所以应尽量增多点数，但也不必过多。适中的DFT点数既可以较好地反映频谱信息，同时又能保证运算效率。【研讨题目3】 基本题3已知一离散序列为 x k=AcosW0k+Bcos ( (W0+DW)k)。用长度N=64的哈明窗对信号截短后近似计算其频谱。试用不同的A和B的取值，确定用哈明窗能分辩的最小的谱峰间隔中c的值。（M2-3） 【题目分析】本题讨论用哈明窗计算序列频谱时的频率分辨率。使用加窗的办法截短时域时会产生频谱混叠，并影响频谱分辨率，对于哈明窗，旁瓣幅度较小可以不考虑频谱泄露。当主瓣宽度等于频谱间隔时，两谱峰一定能分辨开即当，c=4一定能分辨开。而课本中给出的参考值是c=2，因而本题从c=2开始取值取到4，间隔为0.25；(1)先令A=B=1;【仿真程序1与结果】%3A=1;B=1;w0=0.5*pi;N=64;d1=2*2*pi/N;d2=2*2.25*pi/N;d3=2*2.50*pi/N;d4=2*2.75*pi/N;d5=2*3.0*pi/N;d6=2*3.25*pi/N;d7=2*3.5*pi/N;d8=2*3.75*pi/N;d9=2*4.0*pi/N;k=0:N-1;L=16*N;x1=A*cos(w0*k)+B*cos(w0+d1)*k);x2=A*cos(w0*k)+B*cos(w0+d2)*k);x3=A*cos(w0*k)+B*cos(w0+d3)*k);x4=A*cos(w0*k)+B*cos(w0+d4)*k);x5=A*cos(w0*k)+B*cos(w0+d5)*k);x6=A*cos(w0*k)+B*cos(w0+d6)*k);x7=A*cos(w0*k)+B*cos(w0+d7)*k);x8=A*cos(w0*k)+B*cos(w0+d8)*k);x9=A*cos(w0*k)+B