北邮大二下 Matlab实验报告

北京邮电大学

《数字信号处理》Matlab实验

学院：

班级：

姓名：

学号：

班内序号：

实验一

【实验要求】

实现重叠相加和重叠保留算法，完成线性卷积的分段计算。

【实验原理】

一、算法产生背景

DFT是连续傅里叶变换在时域和频域上都离散的形式，将时域信号的采样变换为在离散时间傅里叶变换频域的采样。

对于线性非移变离散系统，可由线性卷积表示时域输入输出关系。即y(n)=x(n)*h(n)。

通常采用循环卷积降低运算量，但实际中往往无法满足对信号处理的实时性要求。因此，产生了重叠相加法和重叠保留法两种典型的算法，用以快速计算线性卷积，成为了DFT的一个重要应用。

二、算法基本思想

1.重叠相加法

重叠相加法是将待过滤的信号分割成长为N的若干段，每一段都可以和有限时宽单位取样响应作卷积，再将过滤后的各段重叠相加。

具体算法实现：建立缓存序列，每次输入N点序列，通过计算x(n)和h(n)的循环卷积实现线性卷积运算，将缓存的M-1点序列和卷积结果相加，并输出前N点作为计算结果，同时缓存后M-1点，如此循环，直至所有分段计算完毕，则输出序列y(n)为最终计算结果。

2.重叠保留法

重叠保留法相当于将x(n)和h(n)作循环卷积，然后找出循环卷积中相当于线性卷积的部分。在这种情况下，将序列y(n)分为长为N的若干段，每个输入段和前一段有M-1个重叠点。此时只需要将发生重叠的前M-1个点舍去，保留重叠的部分并输出，则可获得序列y(n)。

【流程图设计】1.重叠相加法

2.重叠保留法

【MATLAB源代码】

1.重叠相加法

（1）function y = ovrlplus (x,h,N) %重叠相加法实现

M = length(h); %获得h(n)的长度

if N

N = M+1;

end

L = M+N-1; %循环卷积与线性卷积结果相同时需要进行运算的最少点数

Lx = length(x); %获得x(n)的长度

T = floor(Lx/N); %确定分段数T floor 向下取整函数

t = zeros(1,M-1); %初始化序列t(n)

x = [x,zeros(1,(T+1)*N-Lx)]; %不足的分段补零

y = zeros(1,(T+1)*N); %生成输出序列y(n),长度足够长

for i=0:1:T

xi=i*N+1;

x_seg = x(xi:xi+N-1); %选择循环卷积计算时的分段x(n)

y_seg = circonvt(x_seg,h,L); %调用循环卷积计算线性卷积

y_seg(1:M-1) = y_seg(1:M-1)+t(1:M-1);%完成重叠相加

t(1:M-1) = y_seg(N+1:L); %重新对t(n)赋值为保留的后M-1点

y(xi:xi+N-1) = y_seg(1:N); %直接输出前N个点

end

y=y(1:Lx+M-1); %取出最终的输出序列

（2）function[y]=circonvt(x1,x2,N)%循环卷积实现

if length(x1)>N

error('N must be>=the length of x1')

end

if length(x2)>N

error('N must be>=the length of x1')

end

x1=[x1 zeros(1,N-length(x1))];

x2=[x2 zeros(1,N-length(x2))];

m=[0:1:N-1];

x2=x2(mod(-m,N)+1);

H=zeros(N,N);

for n=1:1:N

H(n,:)=cirshftt(x2,n-1,N);

end

y=x1*conj(H');

（3）function y=cirshftt(x,m,N)%循环移位实现

if length(x)>N

error('N must be >= the length of x')

end

x=[x zeros(1,N-length(x))];%补零函数（x(n)的长度小于N，将不够的地方全部补

零）

n=[0:1:N-1]; n=mod(n-m,N);y=x(n+1);

2.重叠保留法

（1）function[y]=ovrlpsav(x,h,N)%实现重叠保留的主函数

Lenx=length(x);

M=length(h);

M1=M-1;

L=N-M1;

h= [h zeros(1,N-M)];

x=[zeros(1,M1),x,zeros(1,N-1)];

K=floor((Lenx+M1-1)/(L));

Y=zeros(K+1,N);

for k=0:K

xk=x(k*L+1:k*L+N);

Y(k+1,:)=circonvt(xk,h,N);

end

Y=Y(:,M:N)';

y=(Y(:))';

（2）function[y]=circonvt(x1,x2,N)%循环卷积

if length(x1)>N

error('N must be>=the length of x1')

end

if length(x2)>N

error('N must be>=the length of x1')

end

x1=[x1 zeros(1,N-length(x1))];

x2=[x2 zeros(1,N-length(x2))];

m=[0:1:N-1];

x2=x2(mod(-m,N)+1);

H=zeros(N,N);

for n=1:1:N

H(n,:)=cirshftt(x2,n-1,N);

end

y=x1*conj(H');

（3）function y=cirshftt(x,m,N)%循环移位

if length(x)>N

error('N must be >= the length of x')

end

x=[x zeros(1,N-length(x))];%补零函数（x(n)的长度小于N，将不够的地方全部补零）

n=[0:1:N-1]; n=mod(n-m,N);y=x(n+1);

（4）function[XK]=dft(xn,N);%DFT函数

n=[0:1:N-1];

k=[0:1:N-1];

WN=exp(-j*2*pi/N);

nk=n'*k;

WNnk=WN.^nk;

XK=xn*WNnk;

【实验结果与分析】1.重叠相加法

x1=[1:10]

x2=[1,0,-1]

ovrlplus (x1,x2,4)

2.重叠保留法

x1=[1:10]

x2=[1,0,-1]

ovrlpsav(x1,x2,4)

【运算量分析】

有限长因果序列x（n）、h（n）的长度分别为Lx和M，直接计算线性卷积y（n），y（n）可视为M个序列的叠加结果，序列长度为Lx，所以每生成一个序列需完成Lx次乘法，共需完成M*Lx次乘法运算。这M个序列依次向右移动一位故需（Lx-1）*（M-1）次加法运算。N越大快此算法的运算量相对来说越小，即实现重叠相加法和重叠保留法的运算量节省的越多。

◎结合教材3.5.1节作运算量分析

1.重叠相加法

分段长度

4 10 100 1000 10000

序列长度

10 0.013151 0.000266 0.018963 0.027378 0.074834

100 0.001173 0.000605 0.000319 0.000917 0.010007 1000 0.009936 0.002824 0.00092 0.000746 0.013927 10000 0.074181 0.027121 0.006682 0.00406 0.009922 100000 0.606933 0.265476 0.052852 0.030106 0.048245 1000000 6.033528 2.483526 0.396206 0.246052 0.33387

10000000 59.93776 25.10466 3.897444 2.328122 3.231499

2.重叠保留法

分段长度

序列长度

4 10 100 1000 10000

10 0.000512 0.000225 0.001249 0.000881 0.008454

100 0.001014 0.000514 0.000272 0.000836 0.008802

1000 0.009297 0.002402 0.000787 0.000687 0.011107 10000 0.05674 0.023194 0.006122 0.003941 0.011083 100000 0.558865 0.216818 0.039175 0.026044 0.042643 1000000 5.221879 2.144089 0.357428 0.236277 0.327539 10000000 52.2464 22.02215 3.593282 2.335725 3.251762 由上表可以看出，同重叠相加法类似，随着数据规模的增大，运算耗时呈线性增长，算法的时间复杂度为O(n)，其中n为数据规模。同样由于分配的缓存空间只由分段长度确定，空间复杂度为O(1)。综合考察，重叠保留法也具有较好的时间和空间复杂度。当数据量达到千万量级时，运算延时最少大约为2.335s，可运用于对信号的实时处理。

【问题与分析】

本次试验中遇到的主要问题就是不熟悉Matlab软件的使用方法和编程语言，经过认真阅读教科书和课程讲义上的Matlab介绍以及例题分析，终于完成代码编写和实验操作。

实验二周期序列的谱分析

【实验要求】

利用DFT分析模拟信号

()()

cos16

a

x t tπ

=

之频谱。

一、设定采样周期T并说明原因；

由奈奎斯特抽样定理可知，应有f0>=2fmax，取f0=20Hz，可保证抽样后的信号不产生混叠。

f0=20Hz,

T0=0.05s,

ω=Ω×T=16π×0.05=0.8π.

二、若令()cos(16)x n nT π=，确定该序列之周期N 并说明原因；

ω=Ω×T=16π×0.05=0.8π, 2π÷ω=2.5,

N 取2.5的倍数10。

三、绘制10个周期内()x n 的取值情况； >>n=[0:99];

>>x=cos(0.8*pi*n); >>plot(n,x)

四、令1()x n 表示()x n 的主值序列，绘制

1(())

DFT x n ，解释取值情况；

>>n=[0:9];

>>x1=cos(0.8*pi*n); >>y=fft(x1,10); >>stem(abs(y));

主值序列可以反映周期序列的全貌。

五、令1()x n 表示()x n 的任一周期，绘制

1(())

DFT x n ，解释取值情况；

>>n=[21:30];

>>x2=cos(0.8*pi*n); >>y=fft(x2); >>stem(abs(y));

任意周期序列的频谱应该和主值序列的频谱相同。

六、令1()x n 表示()x n 的2个周期，绘制

1(())

DFT x n ，解释取值情况;

>>n=[0:19];

>>x1=cos(0.8*pi*n); >>y=fft(x1); >>stem(abs(y));

取两倍周期的DFT ，结果是原点数的两倍，并在其他点补零即可。

七、（选作）若1()()()M x n x n R n =?，而M 不是()x n 周期的整数倍，绘制

1(())

DFT x n ，解释取值情况。

>>n=[14:29];

>>x1=cos(0.8*pi*n); >>y=fft(x1); >>stem(abs(y)); 取M=16，得

1(())

DFT x n 如下图，由图可见出现了失真。

实验三 梳状滤波器的应用

【实验要求】

录制一段自己的话音，时间长度及取样频率自定；对该段声音加入一次反射、三次反射和无穷多次反射。

【实验原理】

回声往往是原始声音衰减后的多个延迟叠加而组成的，因此回声可以用延迟单元来生成。

x(n)表示原始声音信号，α为衰减系数，T 为延迟周期，回声信号

)(......)2()()()(y 2NT n x T n x T n x n x n N -?++-?+-?+=ααα.Z 变换后的系统函数H （Z)可由梳状滤波器实现。

MATLAB filter 函数可用来仿真差分方程，本次实验用的就是这个函数。

【实验过程与结果】

在同一张图上，绘制原声音序列()x n 、加入一次反射后的声音序列1()x n 、加入三次反射后的声音序列3()x n 和加入无穷多次反射后的声音序列()I x n 。 源代码：

>>[x, fs] =wavread(‘C:\Users\Moon\Desktop\mine.wav ’) ; >>sound(x, fs);

>>a = 0.6; T = 0.2;%衰减系数与延时

>>y1 = filter([1, zeros(1,T*fs-1), a], 1, x);%加一次反射 >>sound(y1, fs);

>>y2=filter([1,zeros(1,T*fs-1),a,zeros(1,T*fs-1),a^2],1,x);%加两次反射

>>sound(y2, fs);

>>y3=filter([1,zeros(1,T*fs-1),a,zeros(1,T*fs-1),a^2,zeros(1,T*fs-1), a^3], 1, x); %加三次反射 >>sound(y3, fs);

>>y0 = filter(1, [1, zeros(1,T*fs-1), a], x); %加无数次反射

>>sound(y0, fs);

>>plot(y0, 'y')%加无数次反射波形>>hold on

>>plot(y3, 'r')%加三次反射波形>>hold on

>>plot(y2, 'b')%加两次反射波形>>hold on

>>plot(y1, 'g')%加一次反射波形>>hold on

>>plot(x, 'b')%原音频波形

原音频波形：

加一次反射：

加两次反射：

加三次反射：

加无数次反射：

【效果分析】

结合上述各序列，分析延时、衰减系数对回声效果的影响（提示：定量考察序列

()

x n、

1()

x n、

3()

x n和()

I

x n之间的区别）。

从原始声音信号到无穷次反射信号，信号的幅值按照衰减系数的幂次在衰减，反射次数越高，衰减越厉害。

延时造成信号的水平位移，是回声产生的必要条件。延时只要大于人的听觉0.1s,人能明显感觉到回声。

存在一定的延时，减小衰减系数，回声的效果越好。