信号时域采样-频谱分析(matlab)

基于matlab的时域信号采样及频谱分析一：主要设计方法与步骤：1.画出连续时间信号0sin()()tut-atx(t)=Ae的时域波形及其幅频特性曲线，其中，幅度因子444.128A，衰减因子222.144a，模拟角频率0222.144；2.对信号()xt进行采样，得到采样序列0()sin()()antxnAenTun，050n，其中，1sTf为采样间隔，通过改变采样频率可改变T，画出采样频率分别为200Hz，500Hz，1000Hz时的采样序列波形；3.对不同采样频率下的采样序列进行频谱分析，绘制其幅频和相频曲线，对各频率下采样序列()xn和()xt的幅频曲线有无差别，如有差别说明原因；4.设系统单位抽样响应为5()()hnRn，求解当输入为()xn时的系统响应()yn，画出()xn，()hn，()yn的时域波形及幅频特性曲线，并利用结果验证卷积定理的正确性（此内容将参数设置为444.128A，222.144a，0222.144W，1000fs）；5.用FFT对信号()xn，()hn，()yn进行频谱分析，观察与4中结果有无差别；6.由采样序列()xn恢复出连续时间信号1()xt，画出其时域波形，对比1()xt与原来的连续时间信号()xt的时域波形，计算并记录两者最大误差。二：详细程序及仿真波形分析1.连续时间信号()xt及其200/500/1000HzHzHz频率抽样信号函数()xn%绘制信号x(n)的幅度谱和相位谱clcclearallclosealln=0:50%定义序列的长度是50A=input('请入A的值A:')%设置信号的有关参数a=input('请入a的值a:')w0=input('请入w0的值w0:')T1=0.005T2=0.002T3=0.001T0=0.001x=A*exp(-a*n*T0).*sin(w0*n*T0)y1=A*exp(-a*n*T1).*sin(w0*n*T1)y2=A*exp(-a*n*T2).*sin(w0*n*T2)y3=A*exp(-a*n*T3).*sin(w0*n*T3)closeallsubplot(2,1,1)stem(n,x)%绘制x(n)的图形gridontitle('离散时间信号')subplot(2,1,2)plot(n,x)gridontitle('连续时间信号')05101520253035404550-50050100150离散时间信号05101520253035404550-50050100150连续时间信号figure(2)subplot(3,1,1)stem(n,y1)gridontitle('200Hz理想采样信号序列')subplot(3,1,2)stem(n,y2)gridontitle('500Hz连续时间信号')subplot(3,1,3)stem(n,y3)gridontitle('1000Hz连续时间信号')05101520253035404550-1000100200200Hz理想采样信号序列05101520253035404550-2000200500Hz连续时间信号05101520253035404550-20002001000Hz连续时间信号k=-25:25W=(pi/12.5)*kw=W/piY1=y1*exp(-j*pi/12.5).^(n'*k)figure(3)subplot(2,1,1)plot(w,abs(Y1))gridxlabel('w')ylabel('幅度')title('200Hz理想采样信号序列的幅度谱')axis([-2201000])subplot(2,1,2)plot(w,angle(Y1))gridxlabel('w')ylabel('幅角')title('200Hz理想采样信号序列的相位谱')-2-1.5-1-0.500.511.5205001000w幅度200Hz理想采样信号序列的幅度谱-2.5-2-1.5-1-0.500.511.52-4-2024w幅角200Hz理想采样信号序列的相位谱Y2=y2*(exp(-j*pi/12.5)).^(n'*k)figure(4)subplot(2,1,1)plot(w,abs(Y2))gridxlabel('w')ylabel('幅度')title('500Hz理想采样信号序列的幅度谱')axis([-2201000])subplot(2,1,2)plot(w,angle(Y2))gridxlabel('w')ylabel('幅角')title('500Hz理想采样信号序列的相位谱')-2-1.5-1-0.500.511.5205001000w幅度500Hz理想采样信号序列的幅度谱-2.5-2-1.5-1-0.500.511.52-4-2024w幅角500Hz理想采样信号序列的相位谱Y3=y3*(exp(-j*pi/12.5)).^(n'*k)figure(5)subplot(2,1,1)plot(w,abs(Y3))gridxlabel('w')ylabel('幅度')title('1000Hz理想采样信号序列的幅度谱')axis([-2201000])subplot(2,1,2)plot(w,angle(Y3))gridxlabel('w')ylabel('幅角')title('1000Hz理想采样信号序列的相位谱')-2-1.5-1-0.500.511.5205001000w幅度1000Hz理想采样信号序列的幅度谱-2.5-2-1.5-1-0.500.511.52-4-2024w幅角1000Hz理想采样信号序列的相位谱分析：采样频率为1000Hz时没有失真，500Hz时有横线，产生失真，200Hz时横线加长，失真加大。说明采样频率越大，失真越小。2.设系统单位抽样响应5h()()nRn，求解当输入为()xn时的系统响应()yn，画出()xn，()hn，()yn的时域波形及幅频特性曲线，并利用结果验证卷积定理的正确性（此内容将参数设置为444.128A，222.144a，0222.144W，1000fs）。clcclearallclosealln=1:50%定义序列的长度是50hb=zeros(1,50)%注意：matlab中数组下标从1开始hb(1)=1hb(2)=1hb(3)=1hb(4)=1hb(5)=1closeallsubplot(3,1,1)stem(hb)title('系统hb[n]')m=1:50%设定序列和长度值T=1%设定序列的采样率A=1a=0.4T=1w0=2.0734x=A*exp(-a*m*T).*sin(w0*m*T)subplot(3,1,2)stem(x)title('输入序列x[n]')y=conv(x,hb)subplot(3,1,3)stem(y)title('输出信号y[n]')0510152025303540455000.51系统hb[n]05101520253035404550-101输入序列x[n]0102030405060708090100-101输出信号y[n]figure(2)subplot(3,1,1)plot(n,hb)gridontitle('矩形序列的时域波形')subplot(3,1,2)plot(x)gridontitle('输入信号x[n]的时域波形')subplot(3,1,3)plot(y)gridontitle('输出信号y[n]的时域波形')0510152025303540455000.51矩形序列的时域波形05101520253035404550-101输入信号x[n]的时域波形0102030405060708090100-101输出信号y[n]的时域波形分析：有数字信号处理中经常要进行卷积运算，conv可以用来计算两个有限长序列的卷积，该函数计算的两个序列都是从0n开始。3.用FFT对信号()xn，()hn，()yn进行谱分析，观察与4中结果有无差别。clcclearallclosealln=1:50hb=zeros(1,50)hb(1)=1hb(2)=1hb(3)=1hb(4)=1hb(5)=1closeallsubplot(3,1,1)m=1:50T=1A=1a=0.4T=1w0=2.0734x=A*exp(-a*m*T).*sin(w0*m*T)y=conv(x,hb)subplot(3,1,1)plot(n,abs(fft(hb)))title('h(n)的FFT')subplot(3,1,2)plot(abs(fft(x)))title('x(n)的FFT')subplot(3,1,3)plot(abs(fft(y)))title('y(n)的FFT')051015202530354045500246h(n)的FFT05101520253035404550012x(n)的FFT0102030405060708090100012y(n)的FFT分析：matlab中，计算矢量x的DFT及其逆变换的函数分别为fft和ifft，这两个函数采用了混合算法，当N为质数时，采用的是原始的DFT算法。如果x为一个矩阵时，则调用后计算出每列的N点FFT。4.由采样序列()xn恢复出连续时间信号1()xt，画出其时域波形，对比1()xt与原连续时间信号()xt的时域波形，计算并记录两者最大误差。%设置信号的有关参数clcclearallcloseallA=input('pleaseinputtheA:')a=input('pleaseinputthea:')W0=input('pleaseinputtheW0:')fs=input('pleaseinputthefs:')n=0:49T=1/fst0=10/aDt=1/(5*a)t=0:Dt:t0xa=A*exp(-a*t).*sin(W0*t)K1=50k1=0:1:K1W1max=2*pi*500W1=W1max*k1/K1w1=W1/piXa=xa*exp(-j*t'*W1)x=A*exp(-a*n*T).*sin(W0*n*T)figure(1)subplot(4,1,1)plot(t*1000,xa)title('连续时间信号x(t)')axis([0t0*1000-50150])gridxlabel('t:毫秒')ylabel('x(t)')subplot(4,1,2)plot(w1,abs(Xa))title('连续时间信号频谱Xa(w1)')axis([0100001200])subplot(4,1,3)stem(x)gridxlabel('n')ylabel('x(n)')title('采样序列x(n)')axis([050-15160])x1=spline(n*T,x,t)gridxlabel('t:毫秒')ylabel('x(t)')subplot(4,1,4)plot(t*1000,x1)axis([0t0*10000200])title('由x(n)恢复x1(t)')gridxlabel('t:毫秒')ylabel('x1(t)')axis([045-20160])051015202530354045-50050100150连续时间信号x(t)t:毫秒x(t)0100200300400500600700800900100005001000连续时间信号频谱Xa(w1)05101520253035404550050100150t:毫秒x(t)采样序列x(n)051015202530354045050100150由x(n)恢复x1(t)t:毫秒x1(t)error=max(abs(x1-xa))k2=-25:25W2=(pi/12.5)*k2w2=W2/piX=x*(exp(-j*pi/12.5)).^(n'*k2)%序列的付里叶变换函数figure(2)subplot(2,1,1)plot(w2,abs(X))gridxlabel('w2')yl