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

基于matlab的时域信号采样及频谱分析一：主要设计方法与步骤：画出连续时间信号x(t)=Ae-atsin(t)u(t)0因子A444.128，衰减因子a222.144，模拟角频率 222.144；0x(tx(n)Aeantsin(nT)u(nn5001T 为采样间隔，通过改变采样频率可改变T，画出采样频率分别为200Hz，fs500Hz，1000Hz时的采样序列波形；x(nx(t的幅频曲线有无差别，如有差别说明原因；设系统单位抽样响应为h(n)R5

(n)，求解当输入为x(n)时的系统响应y(n)，画出x(n)，h(n)，y(n)的时域波形及幅频特性曲线，并利用结果验证卷积定理的正确性（此内容将参数设置为A444.128，a222.144，W0222.144，fs1000；用FFT 对信号x(n)，h(n)，y(n)进行频谱分析，观察与4中结果有无差别；x(n)x1

(t)x1

(t)与原来的连续时间信号x(t)的时域波形，计算并记录两者最大误差。二：详细程序及仿真波形分析x(t及其200Hz/500Hz/1000Hzx(n)%绘制信号x(n)的幅度谱和相位谱clcclearallclosealln=0:50%定义序列的长度是50A=input('请入A的值A:')%设置信号的有关参数a=input('请入a的值a:')w0=input('w0T1=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('连续时间信号')15010050

离散时间信号0-500 5 10 15 20 25 30 35 40 45 50连续时间信号150100500-500 5 10 15 20 25 30 35 40 45 50figure(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连续时间信号')200100

理想采样信号序列0-1000 5 10 15 20 25 30 35 40 45 50连续时间信号2000-2000 5 10 15 20 25 30 35 40 45 50连续时间信号2000-2000 5 10 15 20 25 30 35 40 45 50k=-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理想采样信号序列的相位谱')1000度 幅

理想采样信号序列的幅度谱0-2 -1.5

-0.5

0 0.5 w

1.5 2理想采样信号序列的相位谱42角 0幅-2-4-2.5 -2 -1.5 -1 -0.5wY2=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理想采样信号序列的相位谱')

0 0.5 1 1.5 21000度 幅

理想采样信号序列的幅度谱0-2 -1.5

-0.5

0 0.5 w

1.5 2理想采样信号序列的相位谱42角 0幅-2-4-2.5 -2 -1.5 -1 -0.5wY3=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理想采样信号序列的相位谱')

0 0.5 1 1.5 21000度 幅

理想采样信号序列的幅度谱0-2

-1

0 0.5 w

1.5 2理想采样信号序列的相位谱42角 0幅-2-4-2.5 -2 -1.5 -1 -0.5w

0 0.5 1 1.5 21000Hz时没有失真，500Hz时有横线，产生失真，200Hz加长，失真加大。说明采样频率越大，失真越小。设系统单位抽样响应h(n)R(nx(n)y(nx(n)，5h(n)，y(n)的时域波形及幅频特性曲线，并利用结果验证卷积定理的正确性（此内容将参A444.128a222.144，W0222.144fs1000。clcclearallclosen=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('y=conv(x,hb)subplot(3,1,3)stem(y)title('输出信号y[n]')10.5

系统hb[n]00 5 10 15 20 25 30 35 40 45 50输入序列x[n]10-10 5 10 15 20 25 30 35 40 45 50输出信号y[n]10-10 10 20 30 40 50 60 70 80 90 100figure(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]的时域波形')10.5

矩形序列的时域波形00 5 10 15 20 25 30 35 40 45 50输入信号x[n]的时域波形10-10 5 10 15 20 25 30 35 40 45 50输出信号y[n]的时域波形10-10 10 20 30 40 50 60 70 80 90 100分析：有数字信号处理中经常要进行卷积运算，conv可以用来计算两个有限长序列卷积，该函数计算的两个序列都是从n 0开始。用FFT 对信号x(n)，h(n)，y(n)进行谱分析，观察与4中结果有无差别。clcclearallclosen=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')64

h(n)的FFT200 5 10 15 20 25的FFT

30 35 40 45 502100 5 10 15 20 25的FFT

30 35 40 45 502100 10 20 30 40 50 60 70 80 90 100分析：matlabxDFT及其逆变换的函数分别为fft采用了混合算法，当N为质数时，采用的是原始的DFT算法。如果x用后计算出每列的NFFT。x(n)x1

(t)x1

(t)与原连续时间信号x(t)的时域波形，计算并记录两者最大误差。%设置信号的有关参数clcclearallcloseA=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)gridxlabel('t:毫秒')ylabel('x1(t)')axis([045-20160])15010050x 0-500 5

连续时间信号x(t)15 20 25

35 40 4510005000

连续时间信号频谱Xa(w1)150100t 50x 1501001 50x

0 100 200 300 400 500 600采样序列x(n)0 5 10 15 20 25 30由0 5 10 15 20 25

700 800 900 100035 40 45 5030 35 40 45error=max(abs(x1-xa))k2=-25:25W2=(pi/12.5)*k2w2=W2/piX=x*(exp(-j*pi/12.5)).^(n'*kfigu