吉林大学数字信号处理实验报告.doc

数字信号处理课程设计实验报告基础实验实验一 离散时间系统及离散卷积一、 实验目的（1）熟悉MATLAB软件的使用方法。 （2）熟悉系统函数的零极点分布、单位脉冲响应和系统频率响应等概念。 （3）利用MATLAB绘制系统函数的零极点分布图、系统频率响应和单位脉冲响应。 （4）熟悉离散卷积的概念，并利用MATLAB计算离散卷积。 二、实验内容1、离散时间系统的单位脉冲响应（1）选择一个离散时间系统；（2）用笔进行差分方程的递推计算；（3）编制差分方程的递推计算程序；（4）在计算机上实现递推运算；（5）将程序计算结果与笔算的计算结果进行比较，验证程序运行的正确性；2.离散系统的幅频、相频的分析方法（1）给定一个系统的差分方程或单位取样响应；（2）用笔计算几个特殊的幅频、相频的值，画出示意曲线图；（3）编制离散系统的幅频、相频的分析程序；（4）在计算机上进行离散系统的幅频、相频特性，并画出曲线；（5）通过比较，验证程序的正确性；3. 离散卷积的计算（1）选择两个有限长序列，用笔计算其线性卷积；（2）编制有限长序列线性卷积程序；（3）利用计算程序对（1）选择的有限长序列进行卷积运算；（4）比较结果验证程序的正确性。三、实验要求（1）编制实验程序，并给编制程序加注释；（2）按照实验内容项要求完成笔算结果；（3）验证编制程序的正确性，记录实验结果。（4）至少要求一个除参考实例以外的实例，在实验报告中，要描述清楚实例中的系统，并对实验结果进行解释说明。四、实验程序及其结果1、离散时间系统的脉冲响应% y(n)-0.9y(n-1)+ 0.5y(n-2)=x(n)+0.5x(n-1) 设x(n)=2n; 0=n=15clear all;b=1,0.5;a=1,-0.9,0.5;n=0:15; x=impseq(0,0,15)h=filter(b,a,x); % 系统冲击响应subplot(2,1,1);stem(n,h);title(系统冲激响应);xlabel(n);ylabel(h);n=0:15; x=2.*n; nx=0:15;nh=0:15;y=conv_m(x,nx,h,nh);ny=length(y);n=0:ny-1;subplot(2,1,2);stem(n,y);title(系统对x(n)响应);xlabel(n);ylabel(y);B=roots(b);A=roots(a);figure; zplane(B,A);2、离散系统的幅频、相频的分析方法% 差分方程为 y(n)-1.76y(n-1)+1.1829y(n-2)-0.2781y(n-3)% =0.0181x(n)+0.0543x(n-1)+0.0543x(n-2)+0.0181x(n-3)b=0.0181,0.0543,0.0543,0.0181;a=1.000,-1.76,1.1829,-0.2781;m=0:length(b)-1; l=0:length(a)-1; K=500;k=0:1:K;w=pi*k/K; H=(b*exp(-j*m*w)./(a*exp(-j*l*w); %计算频率相应magH=abs(H); % magH为幅度 angH=angle(H); %a ngH为相位subplot(2,1,1); plot(w/pi,magH);grid;xlabel(以pi为单位的频率); ylabel(幅度); title(幅度响应); subplot(2,1,2); plot(w/pi,angH); grid; xlabel(以pi为单位的频率); ylabel(相位); title(相位响应); 3、离散卷积的计算% x=1,4,3,5,3,6,5 , -4=n=2% h=3,2,4,1,5,3, -2=n=3% 求两序列的卷积clear all;x=1,4,3,5,3,6,5; nx=-4:2;h=3,2,4,1,5,3;nh=-2:3;ny=(nx(1)+nh(1):(nx(length(x)+nh(length(h);y=conv(x,h);n=length(ny);subplot(3,1,1);stem(nx,x);xlabel(nx);ylabel(x);subplot(3,1,2);stem(nh,h);xlabel(nh);ylabel(h);subplot(3,1,3);stem(ny,y);xlabel(n);ylabel(x和h的卷积);实验二 离散傅里叶变换与快速傅里叶变换一、实验目的（1）加深理解离散傅立叶变换及快速傅立叶变换概念；（2）学会应用FFT对典型信号进行频谱分析的方法；（3）研究如何利用FFT程序分析确定性时间连续信号；（4）熟悉应用FFT实现两个序列的线性卷积的方法。二、实验内容1、用离散傅立叶变换程序处理时间抽样信号，并根据实序列离散傅立叶变换的对称性，初步判定程序的正确性；2.观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。三角波序列反三角波序列3.已知余弦信号如下 当信号频率 F=50Hz 采样间隔 T=0.000625s ，采样长度 N=64时，对该信号进行傅里叶变换。用FFT程序分析正弦信号，分别在以下情况进行分析。（1）F=50, N=32,T=0.000625；（2）F=50, N=32,T=0.005；（3）F=50, N=32,T=0.0046875；（4）F=50, N=32,T=0.004；（5）F=50, N=64=0.0006254.选定某一时间信号进行N=64点离散傅立叶变换，并且，对同一信号进行快速傅立叶变换，并比较它们的速度。三、实验要求1、编制DFT程序及FFT程序，并比较DFT程序与FFT程序的运行时间。给编制的程序加注释；2.完成实验内容2，并对结果进行分析。在单位圆Z上的变换频谱会相同吗？如果不同，你能说出那个低频分量更多吗？为什么？3.完成实验内容3，并对结果进行分析；4.利用编制的计算卷积的计算程序，给出一下三组函数的卷积结果。 （1） （2）（3） 四、实验程序及其结果分析1、离散傅里叶变换function Xk = dft( xn,N ) %计算N点的DFTn=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnk;End2、（1）三角波 clear all;n=0:1:7;x1=(n+1).*(n=0).*(n=3)+(8-n).*(n=4); % 三角波序列subplot(3,1,1);stem(n,x1); title(三角波序列);xlabel(n);N=8; k=0:N-1;X1=fft(x1,N); magX1=abs(X1);phaX1=angle(X1);subplot(3,1,2);stem(k,magX1);xlabel(k);ylabel(三角波DFT的幅度);subplot(3,1,3);stem(k,phaX1);xlabel(k);ylabel(三角波DFT的相位);（2）反三角波clear all;n=0:1:7;x2=(4-n).*(n=0).*(n=3)+(n-3).*(n=4); % 反三角波序列subplot(3,1,1);stem(n,x2); title(反三角波序列);xlabel(n);N=8; k=0:N-1;X2=fft(x2,N); magX2=abs(X2);phaX2=angle(X2);subplot(3,1,2);stem(k,magX2);xlabel(k);ylabel(幅度);subplot(3,1,3);stem(k,phaX2);xlabel(k);ylabel(相位);结果分析：由图知：三角波和反三角波序列的波形不同，当N=8时，正反三角波的幅频特性相同，因为两者的时域只差一个相位3、余弦信号的FFT分析clea all;F=50; N=32;T=0.000625; % (1) F=50,N=32,T=0.000625n=1:N;x=cos(2*pi*F*n*T);figure(1); subplot(2,1,1); plot(n,x);ylabel(x(n); xlabel(n);title(1) F=50,N=32,T=0.000625);X=fft(x); magX=abs(X);subplot(2,1,2);plot(n,X);ylabel(FFT|X|); xlabel(f(pi);F=50; N=32;T=0.005; % (2) F=50,N=32,T=0.005n=1:N;x=cos(2*pi*F*n*T);figure(2);subplot(2,1,1);plot(n,x);ylabel(x(n);xlabel(n);title(2) F=50,N=32,T=0.005);X=fft(x); magX=abs(X);subplot(2,1,2); plot(n,X); ylabel(FFT|X|); xlabel(f(pi);F=50,N=32,T=0.0046875; % (3) F=50,N=32,T=0.0046875n=1:N;x=cos(2*pi*F*n*T);figure(3);subplot(2,1,1); plot(n,x); ylabel(x(n);xlabel(n);title(3) F=50,N=32,T=0.0046875);X=fft(x);magX=abs(X);subplot(2,1,2);plot(n,X); ylabel(FFT|X|); xlabel(f(pi);F=50,N=32,T=0.004; % (4) F=50,N=32,T=0.004n=1:N; x=cos(2*pi*F*n*T);figure(4);subplot(2,1,1);plot(n,x); ylabel(x(n);xlabel(n);title(4) F=50,N=32,T=0.004);X=fft(x); magX=abs(X);subplot(2,1,2);plot(n,X);ylabel(FFT|X|); xlabel(f(pi); F=50,N=64,T=0.000625; % (5) F=50,N=64,T=0.000625n=1:N;x=cos(2*pi*F*n*T);figure(5); subplot(2,1,1); plot(n,x); ylabel(x(n); xlabel(n);title(5) F=50,N=64,T=0.000625);X=fft(x); magX=abs(X);subplot(2,1,2); plot(n,X); ylabel(FFT|X|); xlabel(f(pi); 结果分析；不同的采样间隔会产生不同的栅栏效应，相当于透过栅栏观赏风景,只能看到频谱的一部分,而其它频率点看不见,因此很可能使一部分有用的频率成分被漏掉，从而产生不同的频谱图，减小栅栏效应可用提高采样间隔也就是频率分辨力的方法来解决。间隔小，频率分辨力高，被“挡住”或丢失的频率成分就会越少。但会增加采样点数，使计算工作量增加。4、卷积计算clear all;n=0:14;x1=1.*(n=0).*(n=14); % 第一组函数 h1=(4/5).n;N=27; y=circonvt(x1,h1,N);disp(y);subplot(3,1,1);i=0:N-1; stem(i,y); ylabel(式（1）卷积);n=0:9;x2=1.*(n=0).*(n=9); % 第二组函数n=0:19;h2=(0.8).*sin(0.5.*n); N=27; y=circonvt(x2,h2,N);disp(y); subplot(3,1,2);i=0:N-1; stem(i,y); ylabel(式（2）卷积);n=0:9; x3=(1-0.1.*n); h3=0.1.*n; % 第三组函数N=17; y=circonvt(x3,h3,N);disp(y); subplot(3,1,3);i=0:N-1; stem(i,y); ylabel(式（3）卷积);实验三 IIR数字滤波器设计一、试验目的1、学习模拟数字变换滤波器的设计方法；2.掌握双线性变换数字滤波器设计方法；3.掌握实现数字滤波器的具体方法。二、实验内容1、设计一个巴特沃思数字低通滤波器，设计指标如下：通带内wwp=0.2pi 幅度衰减不大于1dB；阻带内wws=0.35pi 幅度衰减不小于15dB； 2、编制计算设计的数字滤波器幅度特性和相位特性的程序，并进行实验验证。 3、编制实现该数字滤波器程序并且实现数字滤波 （1）分别让满足所设计的滤波器的通带、过渡带、阻带频率特性的正弦波通过滤波器，验证滤波器性能； （2）改变正弦抽样时间，验证数字低通滤波器的模拟截止频率实抽样时间的函数。三、实验要求1、编制实验内容要求的程序，并给程序加注释；2、根据实验结果，给出自己设计的数字滤波器的幅度特性和相位特性；3.用所设计的滤波器对不同频率的正弦波信号进 行滤波，以说明其特性；4、fp=0.2KHz,Rp=1dB,fs=0.3KHz,As=25dB,T=1ms；分别用脉冲响应不变法及双线性变换法设计一Butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。四、实验程序及结果分析1、巴特沃斯低通滤波器clear all;wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;R=10(-Rp/20);A=10(-As/20);b,a=afd_buttap(wp,ws,Rp,As); C,B,A = sdir2cas(b,a); %将直接型装换成级联型db,mag,pha,w=freqs_m(b,a,0.5*pi); %计算系统频率响应ha,x,t=impulse(b,a);subplot(2,1,1);plot(w/pi,mag);title(幅度响应);grid;ylabel(幅度);xlabel(以pi为单位的频率);subplot(2,1,2);plot(w/pi,pha);title(相位响应);grid;ylabel(相位);xlabel(以pi为单位的频率);2、脉冲响应不变法设计巴特沃斯滤波器% 冲击响应不变法设计数字巴特沃斯低通滤波器clear all;T=0.0001;wp=200*2*pi*T;Rp=1; ws=300*2*pi*T;As=25; %数字滤波器指标 T=1;omegap=wp/T; omegas=ws/T; %转换为模拟域指标cs,ds=afd_buttap(omegap,omegas,Rp,As); %模拟巴特沃斯滤波器的计算% 冲击响应不变法b,a=imp_invr(cs,ds,T); C,B,A=dir2par(b,a); % 将IIR滤波器直接型转换为并联型的函数% 计算数字滤波器频率响应db,mag,pha,grd,w=freqz_m(b,a);subplot(3,1,1);plot(w/pi,mag);title(冲击响应不变法设计巴特沃斯滤波器);axis(0 0.5 0 1);grid; ylabel(|幅度|);xlabel(以pi为单位的频率);subplot(3,1,2);plot(w/pi,db);axis(0 0.8 -400 0);grid;ylabel(对数幅度（db）);xlabel(以pi为单位的频率);3、双线性变换法设计巴特沃斯滤波器clear all;T=0.0001;Fs=1/T;wp=200*2*pi*T;Rp=1; ws=300*2*pi*T;As=25; %数字滤波器指标% 转换成模拟域指标omegap=(2/T)*tan(wp/2);omegas=(2/T)*tan(ws/2);%模拟巴特沃斯滤波器的计算cs,ds=afd_buttap(omegap,omegas,Rp,As);%双线性变换b,a=bilinear(cs,ds,Fs);C,B,A=sdir2cas(b,a);db,mag,pha,grd,w=freqz_m(b,a);subplot(3,1,1);plot(w/pi,mag);ylabel(|幅度|);title(双线性变换法设计巴特沃斯滤波器);xlabel(以pi为单位的频率);axis(0 0.5 0 1);grid;subplot(3,1,2);plot(w/pi,db);axis(0 0.8 -400 0);grid;ylabel(对数幅度（db）);xlabel(以pi为单位的频率);subplot(3,1,3);plot(w/pi,pha);axis(0 0.4 -4 4);grid;ylabel(相位);xlabel(以pi为单位的频率);4、数字滤波clear all; F=0.1;N1=512;T1=0.1;n=0:5:N1;x=sin(2*pi*F*n*T1); wp=0.2*pi;ws=0.35*pi;Rp=1;As=15; T=1;Fs=1/T;OmegaP=(2/T)*tan(wp/2); %预畸变OmegaS=(2/T)*tan(ws/2);ep=sqrt(10(Rp/10)-1); %求平方根Ripple=sqrt(1/(1+ep.2);Attn=1/10(As/20);N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(OmegaP/OmegaS);OmegaC=OmegaP/(10.(Rp/10)-1).(1/(2*N);cs,ds=u_buttap(N,OmegaC);b,a=bilinear(cs,ds,Fs);%括号为模拟滤波器系数，中括号中的为数字滤波器系数mag,db,pha,w=freqz_m(b,a);y=filter(b,a,x);subplot(2,1,1);stem(x);axis(0 100 -1 1);title(输入正弦序列);subplot(2,1,2);stem(y); axis(0 100 -1 1);title(滤波输出序列);结果分析：滤波输出序列前有一段是由于延时造成的。实验四 FIR数字滤波器1、 FIR滤波器的设计clear all;% 用哈明窗函数法来设计FIR滤波器 选N=51wp=0.2*pi; ws=0.35*pi; N=51;n=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(2,2,1);stem(n,hd);title(理想冲激响应);axis(0 N-1 -0.1 0.3);xlabel(n); ylabel(hd(n);subplot(2,2,2);stem(n,h);title(实际冲激响应);axis(0 N-1 -0.1 0.3);xlabel(n); ylabel(h(n);subplot(2,2,3);plot(w/pi,pha);title(相位响应);grid;xlabel(以pi为单位的频率); ylabel(相位);subplot(2,2,4);plot(w/pi,db);title(幅度响应（db）);grid;xlabel(以pi为单位的频率); ylabel(H(db);2、分别用矩形窗、汉宁窗、海明窗、三角窗和 Blackman窗设计一个线性相位的带通滤波器。clear all;w1=0.3*pi; w2=0.5*pi; N=51;n=0:1:N-1;hd=ideal_lp(w2,N)-ideal_lp(w1,N); % 理想低通的冲激响应w_rec=(rectwin(N);h1=hd.*w_rec ; % 矩形窗FIR滤波器的冲激响应db,mag,pha,grd,w=freqz_m(h1,1);figure;plot(w/pi,db);title(矩形窗幅度响应（db）);xlabel(以pi为单位的频率); ylabel(H1(db);w_hann=(hann(N);h2=hd.*w_hann % 汉宁窗FIR滤波器的冲激响应db,mag,pha,grd,w=freqz_m(h2,1);figure;plot(w/pi,db);title(汉宁窗幅度响应（db）);xlabel(以pi为单位的频率); ylabel(H2(db);w_ham=(hamming(N);h3=hd.*w_ham % 哈明窗FIR滤波器的冲激响应db,mag,pha,grd,w=freqz_m(h3,1);figure;plot(w/pi,db);title(海明窗幅度响应（db）);xlabel(以pi为单位的频率); ylabel(H3(db);w_tri=(triang(N);h4=hd.*w_tri % 哈明窗FIR滤波器的冲激响应db,mag,pha,grd,w=freqz_m(h4,1);figure;plot(w/pi,db);title(三角窗幅度响应（db）);xlabel(以pi为单位的频率); ylabel(H4(db);w_bla=(blackman(N);h5=hd.*w_bla % 哈明窗FIR滤波器的冲激响应db,mag,pha,grd,w=freqz_m(h5,1);figure;plot(w/pi,db);title(布莱克曼窗幅度响应（db）);xlabel(以pi为单位的频率); ylabel(H5(db);结果分析：由图可以看出，矩形窗设计的滤波器过渡带最窄，但阻带衰减指标最差，用布莱克曼窗设计的阻带衰减指标最好，但过渡带最宽。 3、FIR滤波器的滤波试验clear all;% 用哈明窗函数法来设计FIR滤波器 选N=51wp=0.2*pi; ws=0.35*pi; N=51;n=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); % 实际的阻带衰减F=0.1; % KHzN1=512;T1=0.1;n=0:5:N1;x=sin(2*pi*F*n*T1);y=filter(h,1,x);subplot(2,1,1);stem(x);title(输入正弦序列);subplot(2,1,2);stem(y);title(滤波输出序列);（1）信号频率低于通带截止频率（2）信号频率位于过渡带频率（2）信号频率大于截止频率结果分析：由图可知，当函数的频率在窗的通带范围内，可以完全滤出来，当频率到达窗的截止频率，则滤出来的波由出现失真。综合实验一、实验目的和要求 （1）实验目的 借助基础实验篇编制的程序，对语音信号进行处理。通过综合实验篇，使得学生能够充分了解信号处理及数字信号处理的过程；使得学生能够对自己设计的频谱分析程序与滤波器程序得到灵活的应用；使得学生能够综合分析信号处理过程各个阶段信号时域与频域特性。 （2）实验要求 1、完成实验内容要求的各项内容。写明设计思路及其设计原理； 2、给出MATLAB代码，并给编制的程序加注释； 3、给出仿真测试结果并对测试结果进行分析；要有仿真验证过程、波形、结果分析、所遇问题及问题解决 。 4、对设计成果做出评价，说明本设计的特点和存在问题，提出改进设计意见。二、实验原理 1、运用MATLAB提供的录音指令进行录音和放音； 2、运用前面实验所做的滤波器对录音信号进行滤波； 3、分别在时域和频域进行对比。三、实验方法与内容 1、录制一段各人自己的语音信号。 2、对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图； 3、给定滤波器的性能指标，采用窗函数法或双线性变换设计滤波器，并画出滤波器的频率响应； 4、用自己设计的滤波器对采集的语音信号进行滤波，画出滤波后信号的时域波形和频谱； 5、对滤波前后的信号进行对比，分析信号的变化； 6、回放语音信号，并与录制语音信号比较； 7、设计出一信号处理系统界面。四、实验程序及结果function pushbutton1_Callback(hObject, eventdata, handles)Fs=12000; %抽样频率time=10; %录音时间x=wavrecord(time*Fs,Fs,double); %录制语音信号wavwrite(x,Fs,D:20101822wanglang.wav); %存储语音信号function pushbutton2_Callback(hObject, eventdata, handles)y=wavread(D:20101822wanglang.wav);wavplay(y,12000);function pushbutton3_Callback(hObject, eventdata, handles) y1,fs,bits=wavread(D:20101822wanglang.wav);y=y1(1:12000); Y=fft(y);axes(handles.axes1); plot(y);axis(5000,1.2*104,-0.2,0.2);title(采样后时域波形);axes(handles.axes2); plot(abs(Y);axis(0,1.5*104,0,200);title(采样后频谱波形);function pushbutton4_Callback(hObject, eventdata, handles)wp=0.2*pi; ws=0.35*pi; N=51;n=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); % 实际的阻带衰减axes(handles.axes3); plot(w/pi,db);title(幅度响应（db）);grid;xlabel(以pi为单位的频率); ylabel(H(db);axis(0 1 -150 20);function pushbutton5_Callback(hObject, eventdata, handles)% 用哈明窗函数法来设计FIR滤波器 选N=51wp=0.2*pi; ws=0.35*pi; N=51;n=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); % 实际的阻带衰减y1,fs,bits=wavread(D:20101822wanglang.wav);y=filter(h,1,y1);Y=fft(y);axes(handles.axes6); plot(y);title(滤波后时域波形);axes(handles.axes5); plot(abs(Y);title(滤波后频谱波形);wavplay(y,7000);实验界面：实验结果及分析提高实验三、基带OFDM系统仿真一 、实验目的借助数字信号处理课程中的快速傅立叶反变换（IFFT）和快速傅立叶变换（FFT）相关知识，实现基带OFDM系统的调制和解调。二、 实验内容1、信号源：首先产生二进制信号，经星座映射为4PSK或16QAM。2、OFDM调制：借助IFFT变换，实现OFDM调制。3、信道模块。OFDM信号通过该模块到达接收端。4、OFDM解调：借助FFT对信道输出信号进行OFDM解调。5、对OFDM解调得到的信号进行星座反映射，还原二进制信号。三、 实验要求1、画出二进制信号、4PSK信号或16QAM信号的星座图。2、画出OFDM符号的时域波形、功率谱图。3、信道为多径+加性高斯白噪声信道。多径信道模型自己选择。4、画出信道输出端信号的时域波形和功率谱图。5、画出OFDM解调后信号的星座图。6、画出整个OFDM系统在不同信噪比下的误码率曲线。四、实验程序及结果分析clear all;carrier_count=200;%子载波数symbols_per_carrier=300; %每子载波含符号数bits_per_symbol=4; %每符号含比特数,16QAM调制IFFT_bin_length=512; %FFT点数PrefixRatio=1/4;%保护间隔与OFDM数据的比例 1/61/4GI=PrefixRatio*IFFT_bin_length ;%每一个OFDM符号添加的循环前缀长度为1/4*IFFT_bin_length 即保护间隔长度为128beta=1/32; %窗函数滚降系数GIP=beta*(IFFT_bin_length+GI); %循环后缀的长度20SNR=15; %信噪比dB%=%=信号产生=baseband_out_length = carrier_count * symbols_per_carrier * bits_per_symbol;%所输入的比特数目carriers = (1:carrier_count) + (floor(IFFT_bin_length/4) - floor(carrier_count/2); %共轭对称子载波映射 复数数据对应的IFFT点坐标conjugate_carriers = IFFT_bin_length - carriers + 2;%共轭对称子载波映射 共轭复数对应的IFFT点坐标baseband_out=round(rand(1,baseband_out_length);%输出待调制的二进制比特流%=16QAM调制=complex_carrier_matrix=qam16(baseband_out);%列向量complex_carrier_matrix=reshape(complex_carrier_matrix,carrier_count,symbols_per_carrier);%symbols_per_carrier*carrier_count 矩阵figure(1);plot(complex_carrier_matrix,*r);%16QAM调制后星座图title(16QAM调制后星座图)axis(-4, 4, -4, 4);grid on%=IFFT=IFFT_modulation=zeros(symbols_per_carrier,IFFT_bin_length);%添0组成IFFT_bin_length IFFT 运算 IFFT_modulation(:,carriers ) = complex_carrier_matrix ;%未添加导频信号 ，子载波映射在此处IFFT_modulation(:,conjugate_carriers ) = conj(complex_carrier_matrix);%共轭复数映射%=signal_after_IFFT=ifft(IFFT_modulation,IFFT_bin_length,2); %OFDM调制 即IFFT变换time_wave_matrix =signal_after_IFFT;%时域波形矩阵，行为每载波所含符号数，列ITTF点数，N个子载波映射在其内，每一行即为一个OFDM符号%=%=添加循环前缀与后缀=XX=zeros(symbols_per_carrier,IFFT_bin_length+GI+GIP);for k=1:symbols_per_carrier; for i=1:IFFT_bin_length; XX(k,i+GI)=signal_after_IFFT(k,i); end for i=1:GI; XX(k,i)=signal_after_IFFT(k,i+IFFT_bin_length-GI);%添加循环前缀 end for j=1:GIP; XX(k,IFFT_bin_length+GI+j)=signal_after_IFFT(k,j);%添加循环后缀 endendtime_wave_matrix_cp=XX;%添加了循环前缀与后缀的时域信号矩阵,此时一个OFDM符号长度为IFFT_bin_length+GI+GIP=660%=OFDM符号加窗=windowed_time_wave_matrix_cp=zeros(1,IFFT_bin_length+GI+GIP);for i = 1:symbols_per_carrier windowed_time_wave_matrix_cp(i,:) = real(time_wave_matrix_cp(i,:).*rcoswindow(beta,IFFT_bin_length+GI);%加窗 升余弦窗end %=生成发送信号，并串变换=windowed_Tx_data=zeros(1,symbols_per_carrier*(IFFT_bin_length+GI)+GIP);windowed_Tx_data(1:IFFT_bin_length+GI+GIP)=windowed_time_wave_matrix_cp(1,:);for i = 1:symbols_per_carrier-1 ; windowed_Tx_data(IFFT_bin_length+GI)*i+1:(IFFT_bin_length+GI)*(i+1)+GIP)=windowed_time_wave_matrix_cp(i+1,:);%并串转换，循环后缀与循环前缀相叠加end%=Tx_data=reshape(windowed_time_wave_matrix_cp,(symbols_per_carrier)*(IFFT_bin_length+GI+GIP),1);%加窗后 循环前缀与后缀不叠加 的串行信号%=加窗的发送信号频谱=symbols_per_average = ceil(symbols_per_carrier/5);%符号数的1/5，10行avg_temp_time = (IFFT_bin_length+GI+GIP)*symbols_per_average;%点数，10行数据，10个符号averages = floor(temp_time1/avg_temp_time);average_fft(1:avg_temp_time) = 0;%分成5段for a = 0:(averages-1) subset_ofdm = Tx_data(a*avg_temp_time)+1):(a+1)*avg_temp_time);%利用循环前缀后缀未叠加的串行加窗信号计算频谱 sub