基于MATLAB的信号分析

1、摘 要本文首先介绍了三种典型数字信号，对离散信号的均值、方差、相关和高斯随机信号的统计特性用MATLAB仿真和分析，用MATLAB实现离散信号的加减运算。其次编程实现了三种典型离散信号的离散傅里叶变换，显示时域信号和频谱图形（幅度值和相位谱），最后用经典功率谱估计中的周期图估计法、Bartlett谱估计法及Welch谱估计法，对正弦序列加高斯随机序列进行功率谱估计，并且用时域提取法进行提取。关键词：离散信号 功率谱估计 提取目 录摘 要1前 言1一 基本原理1二 系统分析3三 详细设计4总 结16附 录17参考文献27致 谢280前 言 数字信号处理无论在理论上还是技术上都有了突破性的发展。它

2、在工业中的成功是由于开发和利用了廉价的硬件和软件，现在各个领域的新工艺和新应用中现在都想利用DSP算法，这就导致了对具有DSP基础的电气工程师的巨大需求，同时也对DSP的电子仿真提出了更高的要求。 随着MATLAB的出现和不断完善，尤其是MATLAB的信号分析工具箱的提出，越来越多的电气工程师门已经意识到用MATLAB来解决电信工程的实际问题是一种省事又省力的选择，由于MATLAB在信号处理与分析等方面有很大的优点，在过去的10年中，MATLAB已经成为DSP应用中分析和设计的主要仿真工具。MATLAB是一个高性能的科学计算软件，广泛应用于数学计算、算法开发、数学建模、系统仿真、数据分析处理及

3、可视化、科学和工程测绘、应用系统开发等。当前它的使用范围涵盖了工业、电子、医学、医疗、建筑等领域。MATLAB用于数字信号处理与分析仿真实验平台的设计已经十分普遍。一 基本原理1.1离散时间信号（离散序列的时域分析）信号可以表示为一个或多个自变量的函数。在信号的分类中，自变量连续变化的信号称为连续时间（模拟）信号，自变量离散变化的信号称为离散时间信号。一个连续（模拟）信号通常用表示其中变量t可代表任何物理量，一般都假定它代表时间，以秒计。一个离散信号通常用表示，其中变量n是整数值并在时间上代表一些离散时刻；因此，它也称作离散时间信号。离散时间信号通常称为(离散)序列。在通用的计算机或专用的数字

4、信号处理器中，只能用有限位二进制数来表示离散时间信号的幅度，即离散时间信号的幅度只能取有限个离散值。这类在幅度上只能取有限个离散值的离散时间信号称为数字信号。在分析离散时间信号时，常将复杂的离散时间信号表示为基本离散序列的线性组合，这些基本离散时间信号也是实际应用中常用的序列。因此基本离散序列在离散时间信号分析中起着非常重要的作用。1.2傅里叶变换（频域分析）离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT的频域采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组

5、序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作其周期延拓的变换。在实际应用中通常采用快速傅里叶变换计算DFT。在各种信号序列中，有限长序列占重要地位。对有限长序列可以利用离散傅立叶变换(DFT)进行分析。DFT不但可以很好的反映序列的频谱特性，而且易于用快速算法(FFT)在计算机上进行分析。有限长序列的DFT是其z变换在单位圆上的等距离采样，或者说是序列傅立叶的等距离采样，因此可以用于序列的谱分析。1.3数字滤波器及滤波滤波器的作用是利用离散时间系统的特性对输入信号波形(或频谱)进行加工处理，或者说利用数字方法按预定的要求对信号进行变换。数字滤波器,

6、是数字信号处理中及其重要的一部分。随着信息时代和数字技术的发展，受到人们越来越多的重视。数字滤波器可以通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊功能。数字滤波器种类很多，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种，即有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。28二 系统分析2.1 系统设计流程：（1） 明确典型数字信号在MATLAB中的产生方法，产生正弦信号、周期信号和随机

7、信号。对确知信号（正弦、周期）进行频域分析，对随机信号进行统计分析（均值、方差、相关）和频域分析。（2） 分别进行正弦信号与随机信号、周期信号与随机信号的加减运算，并对正弦信号与随机信号组成的混合信号进行功率谱分析。系统初始化（3）设计数字低通滤波器对上述产生的混合信号进行滤噪处理，恢复原始正弦信号产生数字信号是否为随机信号？统计分析、频域分析频域分析数字信号的加减运算对正弦余随机组成的混合信号进行功率谱分析设计滤波器对混合信号进行滤波恢复原始正弦信号三 详细设计3.1典型数字信号的产生与分析3.1.1正弦序列 （3-1） 其中是幅度,而是相位（弧度rad）。可用MATLAB函数cos（或si

8、n）产生正弦序列。 下图是正弦信号的正弦序列、正弦序列时域信号及正弦序列幅频、相频特性图，运行MATLAB得到其正弦信号的均值为0.0000 ,方差为2.0328图3-1正弦序列 正弦序列在点DFT，正弦序列的幅频、相频特性图图3-2正弦序列幅频、相频特性3.1.2 周期序列 如果有序列则序列是周期的。满足上面关系的最小整数成为基波周期。现用表示一个周期序列。为了从一个周期的产生的个周期，可将重复次。 下图是周期序列及幅频、相频特性图，其均值为1.5000 方差为1.2500图3-3 周期序列 图3-4周期序列的幅频、相频特性3.1.3 高斯随机序列 利用MATLAB很容易产生两类随机信号：

9、Rand(1,N)在0.1区间上产生N点均匀分布的随机序列； Randn(1,N)产生均值为0，方差为1的高斯随机序列，也就是白噪声序列。 图3-5正弦信号加高斯随机信号自相关3.2 离散信号进行加减运算 若信号： （3-2）值得注意的是，当序列和的长度不等或位置不对应时，首先应使两者的位置对齐，然后通过zeros函数左右补零使其长度相等后再相加。下图是正弦序列加减随机序列图：图3-6离散信号加减运算3.3 经典法功率谱估计 （3-3）经典功率谱估计主要包括Tukey提出的相关函数估计，Schuster提出的周期图估计，间接法（BT法）、Bartlett法谱估计、Welch法谱估计、Butta

10、ll法谱估计，以及Brurg提出的最大熵谱估计。本节主要运用基于直接法的功率谱估计、基于间接法的功率谱估计、Welch法功率谱估计对正弦信号加高斯噪声信号的功率谱估计。3.3.1基于直接法的功率谱估计 直接法是Schuster于1899年提出的，又称周期图法。它将平稳信号序列的N点观察数据视为能量有限的信号，直接取的傅里叶变换，然后再取其幅值的平方，并除以数据长度N，作为序列功率谱的估计。 假设直接估计出的功率谱为，则： 直接法之所以得到广泛使用，是由于它于序列的频谱有对应关系，可以采用FFT算法来快速计算。但是在直接法功率谱估计中，对无限长的平稳信号序列进行截断，这等于对无限长的序列加以矩形

11、窗，使之变成有限长的数据。这也意味着对自相关函数的加窗，使得功率谱与窗函数的卷积。这种频域卷积会产生频谱泄露，容易使弱信号的主瓣被强信号的旁瓣淹没，造成频谱的模糊和失真，使得周期图功率谱的分辨率较低。 在MATLAB中可以利用periodogram实现直接法的功率谱估计，其调用格式为： Pxx,w=periodogram(x) Pxx,w=periodogram(x,window) Pxx,w=periodogram(x,window,nfft,fs)Pxx,=periodogram(x,range)Periodogram()Periodogram函数参数的说明如下：x为进行功率谱估计的输入有

12、限长序列。Window用于指定采用的窗函数，缺省时采用矩形窗，窗长等于输入序列x的长度。对于周期图法，算法未指定窗函数，但有限长序列相当于采用矩形窗。nfft设定FFT算法的长度，一般取2的整次幂，缺省值为256。fs为采样频率，缺省值为1.range用于指定频率的取值范围。若range=twosided,它是复序列的默认选项，其频率取值范围为0,fs,但是如果fs为空向量时，则频率取值范围为0,1,如果没有指定fs,则频率取值范围为0,2;若range=onesided，它是实序列的默认选项，其频率范围为0,fs/2。Pxx为输出的功率谱估计值。f为频率向量。w为归一化的频率向量。 下图用公

13、式法和periodogram实现直接法的功率谱估计曲线，由图可以看出利用公式直接计算的序列的功率谱与利用函数periodogram计算得到的功率谱完全一致。 图3-7正弦序列加高斯随机序列图3-8公式直接计算的功率谱曲线图3-9 periodogram函数计算的功率谱估计曲线3.4 基于改进直接法的功率谱估计 对于功率谱直接法估计，它是通过有限长样本序列的FFT变换来计算平稳信号序列功率谱的，当数据长度N太大时，谱曲线起伏加剧，当数据N长度太小，谱的分辨率又不好。因此必须对周期图法的功率谱估计方法进行改进，这里的改进指的是数据方差的改进。当前主要有两中改进方法：一种是周期图的平滑，即采用间接法

14、估计功率谱；另一种是平均法，它将长度为N的数据分成L段，分别求出每段的功率谱，然后加以平均。下面介绍直接法的改进方法Batlett法、Welch法。3.4.1 Batlett法Batlett平均周期图法的方法是将N点的有限长序列，分成互不重叠的L段数据，每段有M个样本，则有LM=N，即： （3-4）其中，为长度是M的矩形窗。然后求取每一段的M个样本的功率谱估计，即： （3-5）最后求出所有M段数据功率谱的平均值，即： （3-6）并将上式作为整个序列的功率谱估计。随着L的增大，平均周期图的方差趋于零，因此它是功率谱的渐进一致估计。虽然分段平均周期图法功率谱估计可以减小估计误差和波动，但是由于这种

15、方法将长信号分段成短信号，从而使功率谱的分辨率下降。在MATLAB中可以利用函数psd来实现Barlett平均周期图法的功率谱估计，其调用格式为：Pxx=psd(x,nfft,fs,window)Pxx,f=psd(x,nfft,fs,window,noverlap)Pxx,f,F=psd(x,nfft,fs,window,noverlap,p)Psd(x,.,dflag) 在使用psd函数计算Bartlett平均周期图功率谱估计时，还应注意以下几点：(1) 由于Bartlett方法中没有指定窗函数，但对于有限长序列来说，相当于采用矩形窗，所以参数window应设为boxcar。(2) Bar

16、tlett方法中没有指定数据重叠，因而参数noverlap应设为0；(3) 利用psd函数实现Bartlett方法，其分段数L默认为： （3-7） 其中fix表示数值朝零方向取整； (4)参数p为置信概率； (5)参数Pxx为Pxx的p100%置信区间估计。图3-10序列的Bartlett法功率谱估计曲线3.4.2 Welch法Welch法，又称加权交叠平均法，它是对Bartlett法的改进，主要体现在两个方面： 在对序列分段时，允许每段数据有部分交叠； 每段数据可以选择其他的窗函数，不一定要求是矩形窗。在每段数据的功率谱时，Welch法的算法与Bartlett法相同，即： （3-8）其中M为

17、每段数据样本， 它保证由Welch方法得到的功率谱估计是无偏估计。由Welch法得到的平均功率谱为： （3-9）因为Welch法允许各数据段交叠，因而数据段L增加，这样方差得到更大的改善。但是数据的交叠又减小了每一段数据的不相关性，使方差的减小不会达到理论计算的程度。尽管如此，由Welch法估计出来的功率谱是渐进无偏的。由于Welch法在估计序列的功率时，使用了窗函数，由窗函数基本知识可知，采用合适的非矩形窗可以减小信号“频谱泄露”，同时也可以增加谱峰的宽度，从而提高频谱的分辨率。但是在对窗函数进行选择时，有如下的要求： (1) 窗口宽度M远小于样本序列长度N,以排除不可靠的自相关值； (2)

18、 当平稳信号为实过程时，为保证平滑周期图和真实功率谱同样也是偶函数，平滑窗函数必须是实偶对称函数； (3) 平滑窗函数应在m=0处有峰值，并随m绝对值增加而单调下降，使可靠的自相关值有较大的权值； (4) 功率谱是频率的非负函数，由于周期图是非负的，因而要求窗函数的傅里叶变换是非负的。在MATLAB中利用函数psd与pwelch都可实现Welch法的功率谱估计，其方法是一样的，只是部分参数设置有所不同，psd函数的调用格式与Bartlett法功率谱估计相似，这里仅说明函数pwelch的调用格式。Pwelch函数的调用格式如下：Pxx,w=pwelch(x)Pxx,w=pwelch(x,wind

19、ow)Pxx,w=pwelch(x,window,noverlap)Pxx,w=pwelch(x,window,noverlap,nfft)Pxx,f=pwelch(x,window,noverlap,nfft,fs).=pwelch(x,window,noverlap,.range) Pwelch(.) 在使用pwelch函数计算序列功率谱时，还应注意以下两点： (1) 序列的分段数L为： （3-10） (2) 参数range用来指定频率的取值范围。若range=half，频率取值范围为0,Fs/2;若range=whole，取值范围为0,Fs; 在MATLAB中，还可以利用函数csd实现W

20、elch法的互功率谱估计，该函数的调用格式如下： Pxy=csd(x,y) Pxy=csd(x,y,nfft) Pxy,f=csd(x,y,nfft,fs) Pxy=csd(x,y,nfft,fs,window) Pxy=csd(x,y,nfft,fs,window,numoverlap) Pxy=csd(x,y,.dflag) Pxy,Pxyc,f=csd(x,y,nfft,fs,window,numoverlap,p) csd(x,y,.) 图 3-11序列的Welch法功率谱估计曲线3.5 正弦信号的提取将随机信号与正弦信号进行相加，得到混合信号，正弦信号与随机信号如图所示：图3-12正

21、弦信号与随机信号总 结本次设计的核心内容是利用MATLAB，实现三种典型的数字信号，并对实现它们的均值、方差、相关及高斯随机序列高斯随机信号的统计特性的编程，实现连续时间周期信号频域分析的仿真波形，运用经典法功率谱估计对正弦信号加噪进行频谱估计，并用滤波器进行信号的提取。整个设计过程中对信号与系统与数字信号处理有了更深的理解，比如傅立叶变换、信号频谱分析等；其次，实现过程是通过MATLAB软件设计完成的，MATLAB 的图形功能强大，具有良好的人机界面，此次设计过程中熟练了MATLAB的编程，掌握了很多函数的作用及使用方法；最后，通过此次课程设计，我对设计所用到的软件MATLAB有了更加深刻地

22、了解，MATLAB不管在数值计算方面的功能很强大，而且其图形仿真功能更能满足各个领域的需要，因此我们以后更要经常运用MATLAB软件，使其成为自己不可或缺的工具。 附 录%正弦信号k1=-30;k2=30;k=k1:k2;w=pi/6;f=2*sin(k*0.5*w);figure(1);hold on;stem(k,f,'filled');title('正弦序列');xlabel('时间(k)');ylabel('幅值f(k)');fprintf('正弦信号的均值为%.4f 方差为%.4fn',mean(f),v

23、ar(f,1);N=64;k=0:N-1;w=pi/6;f=2*sin(k*0.5*w);f1=fft(f,N);figure(2);hold on;stem(k,f,'filled');xlabel('k');ylabel('幅值f(k)');grid on;title('时域信号');figure(3);hold on;subplot(2,1,1);stem(abs(f1),'filled');ylabel('幅值');title('幅频特性');grid on;subplot(

24、2,1,2);stem(angle(f1),'filled');ylabel('相角');title('相频特性');grid on;%周期信号统计分析k1=0;k2=3;k=k1:k2;Ts=1;f=k*Ts;xtilde=f'*ones(1,8);xtilde=xtilde(:);xtilde=xtilde'subplot(2,1,1);stem(k,f,'filled');title('一个周期');xlabel('时间(k)');ylabel('幅值f(k)'

25、);subplot(2,1,2);stem(xtilde,'filled');title('周期序列');xlabel('时间(k)');ylabel('幅值f(k)');fprintf('周期信号的均值为%.4f 方差为%.4fn',mean(xtilde),var(xtilde,1);k1=0;k2=3;k=k1:k2;Ts=1;f=k*Ts;xtilde=f'*ones(1,8);xtilde=xtilde(:);xtilde=xtilde'ff=fft(xtilde,32);figure(4

26、);hold on;stem(xtilde,'filled');title('周期序列');xlabel('时间(k)');ylabel('幅值f');grid on;figure(5);hold on;subplot(2,1,1);stem(abs(ff),'filled');ylabel('幅值');title('幅频特性');grid on;subplot(2,1,2);stem(angle(ff),'filled');ylabel('相角');

27、title('相频特性');grid on;.随机信号及其分析.k=150;xn=randn(1,k);f_k=fft(xn,k);fprintf('随机信号的均值为%.4f 方差为%.4fn',mean(xn),var(xn,1);figure(1);hold on;stem(xn,'filled');xlabel('时间（k）');ylabel('幅值xn(k)');title('高斯随机信号');grid on;figure(2);hold on;subplot(2,1,1);stem(abs

28、(f_k),'filled');ylabel('幅值');title('幅频特性');grid on;subplot(2,1,2);stem(angle(f_k),'filled');ylabel('相角');title('相频特性');grid on;.离散信号的加减运算.(1)正弦信号+/-随机信号N=48;k=0:N-1;w=pi/10;f=sin(k*w);xn=randn(1,N);x=f+xn;y=f-xn;figure(1);hold on;subplot(2,2,1);stem(k,

29、f,'filled');xlabel('时间（k）');ylabel('幅值f(k)');title('f'); subplot(2,2,2);stem(k,xn,'filled');xlabel('时间（k）');ylabel('幅值xn(k)');title('xn');subplot(2,2,3);stem(k,x,'filled');xlabel('时间（k）');ylabel('幅值x(k)');title(&

30、#39;f+xn');subplot(2,2,4);stem(k,y,'filled');xlabel('时间（k）');ylabel('幅值y(k)');title('f-xn');(2)周期信号+/-随机信号k=30;k1=1:4;Ts=1;f=k1*Ts;xtilde=f'*ones(1,6);xtilde=xtilde(:);xtilde=xtilde'xn=randn(1,k);x=f+xn;y=f-xn;figure(1);hold on;subplot(2,2,1);stem(xtilde,&

31、#39;filled');xlabel('时间（k）');ylabel('幅值f(k)');title('f'); subplot(2,2,2);stem(xn,'filled');xlabel('时间（k）');ylabel('幅值xn(k)');title('xn');subplot(2,2,3);stem(k,x,'filled');xlabel('时间（k）');ylabel('幅值x(k)');title('f

32、+xn');subplot(2,2,4);stem(k,y,'filled');xlabel('时间（k）');ylabel('幅值y(k)');title('f-xn');Fs=1000;%产生有噪信号n=0:1/Fs:1;N=length(n);e=randn(1,N);w0=100*pi;w1=50*pi;xn=sin(2*w0*n)+sin(2*w1*n)+e;%绘制信号波形figure(1)plot(n,abs(xn)xlabel('n')ylabel('xn')title(

33、9;x(n)=sin(2*wo*n)+sin(2*w1*n)+e(n)')%计算序列的DFTnfft=1024;Xk=fft(xn,nfft);%计算序列的PSDPxx1=abs(Xk).2/N;%绘制功率谱图形index=0:round(nfft/2-1);k=index*N/nfft;plot_Pxx1=10*log10(Pxx1(index+1);figure(2)plot(k,plot_Pxx1)xlabel('频率(Hz)')ylabel('功率（dB)')title('公式直接计算的功率谱')%periodogram函数计算的

34、功率谱window=boxcar(length(xn);Pxx2,f=periodogram(xn,window,nfft,Fs);plot_Pxx2=10*log10(Pxx2(index+1);figure(3)plot(k,plot_Pxx2)xlabel('频率(Hz)')ylabel('功率(dB)')title('periodogram函数计算的功率谱')Fs=1000;%产生有噪信号n=0:1/Fs:1;N=length(n);e=randn(1,N);w0=100*pi;w1=50*pi;xn=sin(2*w0*n)+sin(2*

35、w1*n)+e;%绘制信号波形figure(1)plot(n,abs(xn)xlabel('n')ylabel('xn')title('x(n)=sin(2*wo*n)+sin(2*w1*n)+e(n)')%参数设计nfft=1024;window=boxcar(1001);noverlap=0;p=0.9;Pxx,Pxxc=psd(xn,nfft,Fs,window,noverlap,p);%绘制图形index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1);pl

36、ot_Pxxc=10*log10(Pxxc(index+1);%功率谱对数曲线figure(2)plot(k,plot_Pxx)xlabel('频率(Hz)')ylabel('功率（dB)')title('功率谱对数曲线')%置信区间中的功率谱曲线figure(3)plot(k,plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc)xlabel('频率(Hz)')ylabel('功率(dB)')title('置信区间中的功率谱曲线')Fs=1000;%产生有噪

37、信号n=0:1/Fs:1;N=length(n);e=randn(1,N);w0=100*pi;w1=50*pi;xn=sin(2*w0*n)+sin(2*w1*n)+e;%绘制信号波形figure(1)plot(n,abs(xn)xlabel('n')ylabel('xn')title('x(n)=sin(2*wo*n)+sin(2*w1*n)+e(n)')%用Welch平均估计序列的功率谱nfft=1024;window=hamming(100);noverlap=20;range='half'Pxx,f=pwelch(xn,

38、window,noverlap,nfft,Fs,range);plot_Pxx=10*log10(Pxx);figure(2)plot(f,plot_Pxx)xlabel('频率(Hz)')ylabel('功率（dB)')title('Welch法平均功率谱估计')fp=0.05;fs=3.14;rp=0.3;rs=64;wp=2*pi*fp;ws=2*pi*fs; %通带、阻带截止数字频率n,wn=buttord(wp,ws,rp,rs,'s'); %确定巴特沃兹模拟滤波器的阶次z,p,k=buttap(n); %设计归一化巴特沃兹模拟低通滤波器bp,ap