数字信号处理：MATLAB实验代码整理写在前⾯：本⽂中所有的matlab代码已整理，见下载资源“DSP_matlab.zip”。
数字信号处理之MATLAB实验代码整理
1 常⽤的离散时间信号（DSP1_DiscreteTimeSignals.m）
1）单位样本序列
n1 =-2; n2 =5;
n0 =1;
n =[n1:n2];% n1和n2是序列n的上下界
x =[(n-n0)==0];%给处置序列[]，只有当(n-n0)==0时值才为1
stem(n,x);%⾃变量写前⾯
2）单位阶跃序列
n1 =-2; n2 =5;
n0 =1;
n =[n1:n2];
x =[(n-n0)>=0];%只有当(n-n0)>=0时值才为1
stem(n,x);
3）实指数序列
%实现x(n)= a^n
n =[0:10];
x =(0.9).^n;%在含有数组的乘、除、乘⽅运算时需要点乘，表⽰对应运算
stem(n,x);
4）复指数序列
%实现x(n)= e^((a+bj)*n)
a =2;
b =3;% a表⽰实部，b表⽰虚部
n =1:10;
x =exp((a+b*j)*n);
%取实部，幅度谱
figure(1);stem(n,real(x));
%取虚部，幅度谱
figure(2);stem(n,imag(x));
%计算总幅度，画出幅度谱
AmpX =sqrt(real(x).^2+imag(x).^2);%计算总幅度⽅法⼀
AmpX =abs(x);%计算总幅度⽅法⼆
figure(3);stem(n,AmpX);
%计算相位，画出相位谱
AngleX =angle(x);
figure(4);stem(n,AngleX);
5）正弦序列
A1 =3; w1 =0.1*pi; b1 = pi/3;
A2 =2; w2 =0.5*pi; b2 =0;
n =[0:10];
x = A1*cos(w1*n+b1)+ A2*sin(w2*n+b2);
stem(n,x);
6）随机序列
N =10;
x =randn(1,N);%注意中间是逗号，第⼀个参数必须是1，表⽰在[1:N]产⽣标准正态分布随机序列
stem(x);%不⽤输⼊⾃变量n
2 序列运算（DSP2_SequenceOperation.m）
1）信号相加
function[y,n]=sigadd(x1,n1,x2,n2)
%两个序列分别为x1和x2，对应的⾃变量为n1和n2
%运算结果为y，对应的⾃变量为n
n =min(min(n1),min(n2)):max(max(n1),max(n2));%输出⾃变量取输⼊中最⼩和最⼤的
% y1和y2是将输⼊序列x1和x2的⾃变量范围扩⼤到n，因变量不变
y1 =zeros(1,length(n));%1×length(n)的矩阵，即⼀维数组
y2 = y1;
y1(find((n>=min(n1))&(n<=max(n1))))=x1;%出所有n中对应n1的⾃变量（返回⼀个数组），将因变量x1赋值给y1 y2(find((n>=min(n2))&(n<=max(n2))))=x2;%出所有n中对应n2的⾃变量（返回⼀个数组），将因变量x2赋值给y2 y = y1 + y2;%对于⾃变量相同的两个序列，可以直接相加
end
2）信号相乘
function [y,n]=sigmul(x1,n1,x2,n2)
%两个序列分别为x1和x2，对应的⾃变量为n1和n2
%运算结果为y，对应的⾃变量为n
n =min(min(n1),min(n2)):max(max(n1),max(n2));
y1 =zeros(1,length(n));
y2 = y1;
y1(find((n>=min(n1))&(n<=max(n1))))=x1;
y2(find((n>=min(n2))&(n<=max(n2))))=x2;
y = y1.*y2;%对于⾃闭那辆相同的两个序列，可以直接相乘，注意矩阵和矩阵的对应相乘，要⽤点乘.
end
3）信号加权
x1 =[1,2,3,4,5,6];
a =3;
x2 = a*x1;%直接乘以系数，如果系数是常数不⽤点乘
stem(x2);
4）信号移位
n =[1:10];
x1 = n;
stem(n,x1(n));title('x(n)')
n0 =2;%移位的⼤⼩
n2 = n-n0;%相当于把⾃变量左移n0，由于作图时是把⾃变量和因变量数组对齐的，所以⾃变量左移相当于把信号左移stem(n,x2);title('x(n-n0)')%整个信号左移，移位的结果是x(n+n0)：左加右减
5）信号反转
function[y,n]=sigfold(x,n)
%输⼊序列为x，⾃变量为n；输出序列为y，⾃变量为n
%注意反褶的结果只是对n进⾏反褶：若原信号是x(n-5)，反褶的结果是x(-n-5)，不是x(-n+5) y =fliplr(x);% fliplr将矩阵沿y轴反褶
n =-fliplr(n);%注意反转时n要取相反数
end
6）信号卷积
function [ y,ny ]=conv_m( x,nx,h,nh )
%计算两个离散信号x和h的卷积，对应的⾃变量是nx和nh
%输出卷积结果为y，对应的⾃变量为ny
nyb =nx(1)+nh(1);%卷积结果的⾃变量最⼩值为输⼊序列⾃变量最⼩值之和
nye =nx(length(x))+nh(length(h));%卷积结果地⾃变量最⼤值为输⼊序列⾃变量最⼤值之和ny =[nyb:nye];%确定上下界，得到卷积结果的⾃变量
y =conv(x,h);%⽤conv函数直接计算x和h的卷积（不看⾃变量）
end
7）实验⼀_1：对信号进⾏加权、平移、反褶等基本运算操作
x =[1,2,3,4,5,6,7,6,5,4,3,2,1];%原信号序列
n =[-1:11];
%%画出原信号
subplot(2,2,1)
stem(n,x,'r-');title('x(n)');%⾸先画出原信号x(n)
%%计算0.2*x(5-n)matlab傅里叶变换的幅度谱和相位谱
x1 =0.2*x;%⾸先对信号进⾏加权，得到0.2*x(n)
n1 = n -5;%然后对信号进⾏平移，要得到0.2*x(n+5)，信号应该左移5，所以⾃变量-5
x1 =fliplr(x1);%最后对信号进⾏反褶。注意反褶只对⾃变量操作，得到0.2*x(-n+5)
n1 =-fliplr(n1);
subplot(2,2,2)
stem(n1,x1,'b-');title('0.2*x(5-n)');
%%计算0.3*x(n-3)
x2 =0.3*x;%先加权，得到0.2*x(n)
n2 = n+3;%然后平移，要得到0.3*x(n-3)，左加右减，信号应该右移3，所以⾃变量要+3 subplot(2,2,3)
stem(n2,x2,'g-');title('0.3*x(n-3)');
%%计算y(n)=0.2x(5-n)+0.3x(n)x(n-3)
% y =0.2.*x1 +0.3.*x*x2;
[y1,m1]=sigmul(x,n,x2,n2);%先得到0.3x(n)x(n-3)
[y,n]=sigadd(x1,n1,y1,m1);%再把两个加起来，得到y(n)
subplot(2,2,4)
stem(n,y);title('y(n)');
8）实验⼀_2（主要看⾃编卷积的过程）：
卷积原理：
%%⾃⼰求卷积过程
%先求卷积结果ny
nx =n(1):n(length(n));% x的⾃变量向量nx
nh =n(1):n(length(n));% h的⾃变量向量nh
nyb =nx(1)+nh(1);
nye =nx(length(nx))+nh(length(nh));
ny =[nyb:nye];%卷积结果y的⾃变量向量ny
%再求卷积结果的因变量取值y
%根据卷积原理，这⾥⽤k代替公式中的输出y的⾃变量n,⽤m代表输⼊x的⾃变量m，⽤t代表输⼊h的⾃变量n-m
%其中index_k,index_m,index_t分别是⾃变量向量k,m,t的下标索引
for k =ny(1):ny(length(ny))% k是输出信号y的⾃变量ny的值遍历
index_k = k - nyb +1;%值k对应⾃变量向量ny中的下标索引
y(index_k)=0;%之后要累加，这⾥先初始化为0
for m =nx(1):nx(length(nx))% m是输⼊信号x的⾃变量nx的值的遍历
index_m = m -nx(1)+1;%值m对应向量nx中的下标索引
t = k - m;%对应公式中的n-m，是输⼊信号h的⾃变量nh的值的遍历
if t>=nh(1)&& t<=nh(length(nh))%判断t = k-m是否越界（是否在h的⾃变量的值nh范围内）
index_t = t -nh(1)+1;%值t对应向量nh的下标索引
y(index_k)=x(index_m)*h(index_t)+y(index_k);%通过下标索引到x(m)*h(t)，⽤迭代实现累加
end
end
end
y  %求出的卷积结果是y，对应的⾃变量是ny
3 连续时间信号（见DSP3_ContinuousTimeSignalsSpectrumAnalysis.m）
1）周期信号傅⾥叶级数（以周期⽅波信号的求解为例）
%信号在时域上是周期连续的，傅⾥叶级数在频域上是离散周期的
%先画出周期⽅波信号
T =2;%周期为T
dt =0.00001;
t =-2:dt:2;%为了模拟时间上的连续，引⼊很⼩的变化值dt（如果没有dt，t只有5个值，是离散的）
x =u(t)-u(t-1-dt);%在0和1之间是⾼电平，x只是⼀个周期内的⽅波信号，为后⾯的傅⾥叶级数求解做准
备。注意这⾥要-dt！！不然相位谱会有负值subplot(2,2,1);plot(t,x);%画出⼀个周期内的效果x，连续信号作图⽤plot
x1 =0;% x1是周期性的⽅波。注意要记得赋初值0
for m =-1:1%周期延拓
x1 = x1 +u(t-m*T)-u(t-1-m*T-dt);%在m*T和m*T+1是⾼电平。注意要把x1加上去，不然是覆盖。同样要-dt
end
subplot(2,2,2);plot(t,x1)%画出多个周期的效果x1
%求周期信号的傅⾥叶级数
w0 =2*pi/T;%根据周期计算⾓频率
N =10;%谐波次数
for k =-N:N  %画出正负对称频谱，k即频域上频点的取值
ak(N+1+k)=(1/T)*x*exp(-j*k*w0*t')*dt;
% t'相当于进⾏了求和运算：x是1×length(x)的矩阵，t'是length(x)×1的矩阵
% dt相当于把积分区域分成了许多个⼩的长⽅形进⾏离散求解
%注意这⾥x是⼀个周期内的（傅⾥叶级数的计算公式中是对⼀个周期求解），不能⽤x1
% ak是谐波的权重，当谐波个数N⾜够⼤时，就可以近似地表⽰x(t)
end
subplot(2,2,3);stem([-N:N],abs(ak));%画出幅度谱，离散频谱⽤stem
phi =angle(ak);%计算相位phi
subplot(2,2,4);stem([-N:N],phi);%画出相位谱，离散频谱⽤stem
2）⾮周期信号的傅⾥叶级数（以⼀个[-2,2]上的分段函数为例）
%⾮周期连续信号的傅⾥叶变换（频域上）是连续周期的
dt =0.01;%为了模拟连续时间，引⼊dt
dw =0.1;%为了模拟连续频谱，引⼊dw
t =-10:dt:10;
w =-4*pi:dw:4*pi;
%分段函数
y1 = t+2;%当-2<=t & t<-1时
y2 =1;%当-1<=t & t<1时
y3 =-t+2;%当1<=t & t<2时
y = y1.*(-2<=t&t<-1)+ y2.*(-1<=t&t<1)+ y3.*(1<=t&t<2);%分段函数的实现：⼩括号内为处置序列，true返回1，false返回0 x = y.*(u(t+2)-u(t-2));% x为有限长的信号，只取-2到2
subplot(2,2,1.5);
plot(t,x);title('时域：x(t)'),axis([-5,5,-1,1]),xlabel('t')
grid on
%for k =-4*pi:dw:4*pi
%X(floor((k+4*pi)/dw+1))= x*exp(-j*t'*k)*dt;%类⽐上⾯的傅⾥叶级数
% end
X = x*exp(-j*t'*w)*dt;%第⼆种实现⽅法（求傅⾥叶级数也可以这么实现）
%先通过t'*w实现⼀个length(x)×length(x)的矩阵，和x相乘后得到1×length(x)的矩阵，即傅⾥叶变换结果X
%注意这⾥的t'*w不能写成w*t'（矩阵乘法不满⾜交换律）
subplot(223)
plot(w,abs(X));xlabel('w');title('幅度谱')
grid on
subplot(224)
plot(w,angle(X));xlabel('w');title('相位谱')
3）实验⼆_1：
计算傅⾥叶级数后，利⽤傅⾥叶级数合成信号。并观察谐波次数N对合成信号的影响
T =2;
dt =0.00001;
t =-2:dt:2;
x0 =u(t)-u(t-1-dt);%⼀个周期的信号。⼀定要-dt
%%画出原信号
x =0;
%周期延拓
for m =-1:1
x = x +u(t-m*T)-u(t-m*T-1-dt);
end
subplot(221);plot(t,x);
N =10;%谐波次数
w0 =2*pi/T;
%%
for k =-N:N
ak(N+1+k)=(1/T)*x0*exp(-j*k*w0*t')*dt;%注意要乘以w0⾓频率
end
subplot(222);stem([-N:N],abs(ak));
%%画幅度谱
phi =angle(ak);
subplot(223);stem([-N:N],phi);
%%画相位谱
x2 =zeros(1,length(t));
for m =-2:dt:2
x2(round((m+2)/dt)+1)=0;%注意这⾥要⽤round四舍五⼊！不能⽤floor向下取整
for k =-N:N  %求每⼀个t对应的x2(t)：将各次谐波[-N:N]的信号加起来
x2(round((m+2)/dt)+1)=x2(round((m+2)/dt)+1)+ak(k+N+1)*exp(j*k*w0*m);%通过合成信号的公式合成
end
end
subplot(224);plot(t,x2);%画出合成信号。谐波次数N越⼤，效果越好
4 Z变换（见DSP4_ZTransform.m）

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。