数字信号处理实验报告MATLAB

1、实验内容

一、编写子函数

计算长度为N的序列x(n) （ 0≤n ≤ N-1）的离散时间傅里叶变换，将频率均匀离散化， 一个周期内有M个点。要求画出虚部、实部、幅度、相位，并标注坐标轴。

二、对矩形序列x(n)=RN(n)

1. 用公式表示x(n)的频谱，求出其幅度谱和相位谱； 2. 利用编写的子函数，计算并画出x(n)的频谱

1）固定M，改变N，观察N的取值对频谱的最大值、过零点、第一旁瓣幅度与最大值的比值以及相位谱的影响；

2）固定N，改变M，观察M的取值对幅度谱和相位谱的影响。 如： M=4，26，100 N=4，26，100

三、利用子函数，画出信号x(n)=sin(pi*n/5)和 x(n)=cos(pi*n/4)+cos(pi*n/8) （ 0≤n ≤ N-1）的幅度谱和相位谱。 N分别取为8，16，20，64，75，128，M=256。观察N取不同值时信号频谱的相同和不同之处，为什么会有这样的结果。

2、编程原理、思路和公式

（1）、给定长度为N的序列x(n)的离散时间傅里叶变换（DTFT）的公式为

X(e)x(n)ejwnjwn0N1

x(n)12wjwnX(je)edw。可以看出x(n)的DTFT仍然

是一个连续函数，所以需要将数字角频域w离散化，设一个频率周期内离散点有M个，则第k个点所代表的数字角频率

X(e)X(ejwj2kMw2kM。这样x(n)的DTFT变成：

)x(n)en0N1j2nkM 0kM1，

（2）、因为0≤.n ≤.N-1, 0kM1 设置两层for循环，用sum表示累加求和

X(e)X(ejwj2kM的值，即可实现

)x(n)en0N1j2nkM

（3）对矩形序列x(n)=R4(n)，首先固定M不变，改变N的取值。主函数x=ones(1,N)；M=100; w=0:2*pi/M:2*pi*(M-1)/M;然后调用已经写好

的dtft子函数，求出相位谱和幅度谱，然后调用plot函数即可画出图形。 （4）再次，固定N不变，改变M的取值。x=ones(1,N)；M=4;

w=0:2*pi/M:2*pi*(M-1)/M; 然后调用已经写好的dtft子函数，求出相位谱和幅度谱，然后调用plot函数即可画出图形。

（5）实验三 主函数M=256;N=8;n=0:N-1;x=sin(pi*n/5);

x=cos(pi*n/4)+cos(pi*n/8);然后调用子函数dtft,求出相位谱和幅度谱，然后调用plot函数即可画出图形。

(6)改变N的取值，其余不变，即可观察信号频谱的不同之处。

3、程序脚本，并注释 子函数：

function [Xm,Xp]=dtft(x,M) N=length(x); Xk=zeros(1,M); for k=0:M-1 w=2*pi*k/M; sum=0;

for n=0:N-1

sum=sum+x(n+1)*exp(-j*w*n); end

Xk(k+1)=sum; end

Xm=abs(Xk); Xp=phase(Xk);

实验1

function DTFT(xn,N,M); N=10; n=0:N-1; xn=0:N-1; M=100; k=0:M-1; w=2*pi/M*k;

Xw=xn*exp(-j*(n'*w)); Xw_real=real(Xw);

Xw_imag=imag(Xw); Xw_abs=abs(Xw);

Xw_angle=angle(Xw); subplot(2,2,1); plot(w,Xw_real); xlabel('w'); ylabel('real');

title('DTFT-real'); subplot(2,2,2); plot(w,Xw_abs); xlabel('w'); ylabel('abs');

title('DTFT-abs'); subplot(2,2,3); plot(w,Xw_angle); xlabel('w');

ylabel('angle');

title('DTFT-angle'); subplot(2,2,4); plot(w,Xw_imag); xlabel('w'); ylabel('imag');

title('DTFT-imag');

实验2

x=ones(1,4);M=100;w=0:2*pi/M:2*pi*(M-1)/M; [Xm,Xp]=dtft(x,M);

subplot(2,5,1);plot(w,Xm);

xlabel('w');ylabel('Xm');title('幅频特性');grid; subplot(2,5,6);plot(w,Xp);

xlabel('w');ylabel('Xp');title('相频特性');grid; x=ones(1,26);M=100;w=0:2*pi/M:2*pi*(M-1)/M; [Xm,Xp]=dtft(x,M);

subplot(2,5,2);plot(w,Xm);

xlabel('w');ylabel('Xm');title('幅频特性');grid; subplot(2,5,7);plot(w,Xp);

xlabel('w');ylabel('Xp');title('相频特性');grid; x=ones(1,100);M=100;w=0:2*pi/M:2*pi*(M-1)/M; [Xm,Xp]=dtft(x,M);

subplot(2,5,3);plot(w,Xm);

xlabel('w');ylabel('Xm');title('幅频特性');grid; subplot(2,5,8);plot(w,Xp);

xlabel('w');ylabel('Xp');title('相频特性');grid; x=ones(1,100);M=26;w=0:2*pi/M:2*pi*(M-1)/M; [Xm,Xp]=dtft(x,M);

subplot(2,5,4);plot(w,Xm);

xlabel('w');ylabel('Xm');title('幅频特性');grid; subplot(2,5,9);plot(w,Xp);

xlabel('w');ylabel('Xp');title('相频特性');grid; x=ones(1,100);M=4;w=0:2*pi/M:2*pi*(M-1)/M; [Xm,Xp]=dtft(x,M);

subplot(2,5,5);plot(w,Xm);

xlabel('w');ylabel('Xm');title('幅频特性');grid; subplot(2,5,10);plot(w,Xp);

xlabel('w');ylabel('Xp');title('相频特性');grid;

实验3

M=256;N=8;n=0:N-1; x=sin(pi*n/5);

[Xm,Xp]=dtft(x,M);

subplot(2,2,1);plot(w,Xm);xlabel('w');ylabel('Xm');title('幅频特性');grid;

subplot(2,2,3);plot(w,Xp);xlabel('w');ylabel('Xp');title('相频特性');grid;

x=cos(pi*n/4)+cos(pi*n/8); [Xm,Xp]=dtft(x,M);

subplot(2,2,2);plot(w,Xm);xlabel('w');ylabel('Xm');title('幅频特性');grid;

subplot(2,2,4);plot(w,Xp);xlabel('w');ylabel('Xp');title('相频特性');grid;

4、仿真结果、图形 实验1

实验2

实验3 N=8

N=16

N=20

N=64

N=75

N=128

5、结果分析和结论

（1） 对给定长度为N的序列x(n)的离散时间傅里叶变换（DTFT）将频率均匀离

散化， 一个周期内有M个点，画图可得其幅频特性相频特性。 （2） 根据x(n)的频谱,可以得到：固定M，改变N，当N增加时频谱的最大值

在增加；过零点的数目增多；第一旁瓣幅度与最大值的比值不变。相位谱的变化更加迅速并且逐渐变得杂乱。当M=N时，抽得都是0值附近的数，实部和虚部很小。

（3） 固定N，改变M，可以看出幅度谱的峰值没有太大变化，随着M的增加

旁瓣的数目增多，主瓣的宽度变窄，幅度谱和相位谱整体趋向平稳。 （4） 信号x(n)=sin(pi*n/5)和 x(n)=cos(pi*n/4)+cos(pi*n/8) （ 0≤n ≤ N-1）的幅度

谱和相位谱，N的取值越大信号频谱的最大值越大，尖峰宽度越窄，尖峰越陡峭，边缘振荡越小。

5、遇到的问题、解决方法及收获

（1） 不知道如何写主函数Xk=zeros(1,M); zeros(1,M)表示1行M列个0 （2） 虽然会写离散傅里叶变换公式，但是写两个FOR循环难度较大，写子函数

比较困难。

（3） sum=sum+x(n+1)*exp(-j*w*n);中的x(n+1)是因为MATLAB中，n

的取值是从1开始。

（4） 调用子函数时，应当将该程序和子函数程序放在一个文件夹里，子函数程

序的文件名应该是dtft。