Matlab小波工具箱的使用

-Matlab小波工具箱的使用1 (2011-11-10 20:12:39)转载标签：小波分析连续小波变换尺度系数信号分类：学科知识最近想尝试一下小波的用法，就这matlab的帮助尝试了一下它的例子，顺便翻译了一下帮助的内容，发现matlab帮助做的确实不错，浅显易懂！现把翻译的文档写出来吧，想学习的共同学习吧！小波工具箱简介小波工具箱包含了图像化的工具和命令行函数，它可以实现如下功能：l 测试、探索小波和小波包的特性l 测试信号的统计特性和信号的组分l 对一维信号执行连续小波变换l 对一维、二维信号执行离散小波分析和综合l 对一维、二维信号执行小波包分解（参见帮助Using Wavelet Packets）l 对信号或图像进行压缩、去噪另外，工具箱使用户更方便的展示数据。用户可以做如下选择：l 显示哪个信号l 放大感兴趣的区域l 配色设计来显示小波系数细节工具箱可以方便的导入、导出信息到磁盘或matlab工作空间。具体详见File Menu Options一维连续小波分析这一部分来测试连续小波分析的特性。连续小波分析只需要一个小波函数cwt。在这一部分将学到如下内容：l 加载信号l 对信号执行连续小波变换l 绘制小波系数l 绘制指定尺度的小波系数l 绘制整个尺度小波系数中的最大值l 选择显示方式l 在尺度和伪频率之间切换l 细节放大l 在普通或绝对模式下显示系数l 选择执行小波分析的尺度使用命令行执行连续小波分析这个例子是一个包含噪声的正弦波1. 加载信号load noissin可以使用whos显示信号信息whosNameSizeBytesClassnoissin1x10008000double2. 执行连续小波变换c = cwt(noissin,1:48,db4);函数cwt的参数分别为分析的信号、分析的尺度和使用的小波。返回值c包含了在各尺度下的小波系数。对于这里，c是一个48x1000的矩阵，每一行与一个尺度相关。3. 绘制小波系数cwt函数可以接受第四个参数，来指定函数在执行结束后是否绘制连续小波变换系数的绝对值。另外还可以接受更多的参数来定义显示的不同特性，详见cwt函数。如下面的语句绘制系数结果c = cwt(noissin,1:48,db4,plot);4. 选择分析的尺度cwt函数的第二个参数可以设定任意小波分析的尺度，只要这些尺度满足如下要求l 所有尺幅必须为正实数l 尺度的增量必须为正l 最高的尺度不能超过由信号决定的一个最大值如下面的代码可以执行从2开始的偶数尺度计算c = cwt(noissin,2:2:128,db4,plot);显示结果如下这幅图像很明确的表示出了信号的周期性。使用图形接口做连续小波分析1. 开启一维连续小波工具，只需输入如下命令wavemenu出现如下小波工具箱主菜单选择Continuous Wavelet 1-D菜单项，出现如下一维信号分析连续小波分析工具2. 加载信号选择菜单File-Load Signal，在Load Signal对话框里选择noissin.mat文件，它在matlab安装目录的toolbox/wavelet/wavedemo文件夹下，点击OK加载信号。一维连续小波工具开始加载信号，加载后默认采样频率为1s。3. 执行连续小波变换下面来测试使用db4小波对尺度1到48做小波分析，设置如下4. 点击Analyze按钮在短暂的计算后，工具将绘制小波系数，并在Coefficients line坐标系中绘制尺度为24的小波系数，在local maxima坐标系中绘制各尺度的小波系数最大值。5. 查看小波Coefficients Line在小波系数图中右键点击可以选择展示其他尺度的小波系数，选择后点击New Coefficients Line按钮，Coefficients Line会相应更新。6. 查看Maxima Line点击Refresh Maxima Line按钮，可以显示从尺度1到所选尺度的小波系数的最大值。注意当在系数图中按下鼠标右键并移动时，会在最下面的Info框中显示当前鼠标位于的X位置和尺度。7. 在尺度和伪频率之间切换在右边选择Frequencies，当再在系数图中选择时，在Info中显示的将是Hz。而关于尺度和频率的转换关系，可以看How to Connect Scale to Frequency?8. 选择要显示的坐标系9. 放大细节在系数框中按鼠标左键可以选择放大的范围。10. 选择好放大范围后点击最下面的按钮可以实现指定的放大11. 显示普通系数或系数绝对值两种显示方式的区别在于，普通模式下，颜色映射是在系数的最大最小之间；而绝对模式，颜色映射是在0和最大的系数绝对值之间。图形接口的导入导出信息导入信号到一维连续小波工具首先将要处理的信号保存到mat文件中，要求信号是一维的向量。然后使用工具的File-Load Signal菜单功能，选择此信号文件即可导入信号。文件中第一个一维变量被认为是信号，变量在文件中顺序是按字母排序的。保存小波系数小波分析完成后，点击File-Save-Coefficients，可以将分析结果保存到mat文件。保存后，可以使用load函数加载数据，会看到保存的变量有小波系数coeff、尺度scales、小波的名字wname。一维复信号连续小波分析对于复信号连续小波分析，小波工具箱中对应的函数是cwt。使用命令行分析复信号这个例子是一个带尖头的信号，如下1. 加载信号load cuspamax文件中包含两个变量，caption和cuspamax，前者是此信号的定义，如下caption =x = linspace(0,1,1024); y = exp(-128*(x-0.3).2)-3*(abs(x-0.7).0.4);2. 执行连续小波变换c = cwt(cuspamax,1:2:64,cgau4);3. 绘制结果c = cwt(cuspamax,1:2:64,cgau4,plot);结果显示了四副图形，分别为系数的实部、虚部、模、相角。使用图形接口分析复信号与实信号不同的是，选择Complex Continuous Wavelet 1-D，得到的结果如下具体操作过程与实信号的相似，如下Matlab小波工具箱的使用2 (2011-11-11 09:32:57)转载标签：信号阈值分解系数小波分析分类：学科知识一维离散小波分析工具箱提供了如下函数做一维信号分析：Function NamePurpose分解函数dwt一层分解wavedec分解wmaxlev最大小波分解层数重构函数idwt一层重构waverec全重构wrcoef有选择性重构upcoef单一重构分解结构工具detcoef细节系数抽取appcoef近似系数抽取upwlev分解结构重排去噪和压缩ddencmp使用默认进行去噪和压缩wbmpen使用处罚门槛进行一维或二维去噪wdcbm使用处罚门槛进行一维（使用Birg-Massart方法）wdencmp小波去噪和压缩wden自适应小波去噪wthrmngr门槛设置管理在这一部分，可以学到l 加载信号l 执行一层小波分解l 从系数重建近似和细节l 显示近似和细节l 通过逆小波变换重建信号l 执行多层小波分解l 抽取近似系数和细节系数l 重构第三层近似l 重构第1、2、3层细节l 显示多层分解的结果l 从第三层分解重构原始信号l 从信号中去除噪声l 改善分析l 压缩信号l 显示信号的统计信息和直方图一维分析-使用命令行这个例子包含一个真实世界的信号-测量3天的电功耗。这个信号很典型，因为它包含一个明显的测量噪声，而小波分析可以有效的移除噪声。1. 加载信号load leleccum截取信号s = leleccum(1:3920);l_s = length(s);2. 对信号执行一层小波分解使用db1小波执行一层小波分解，执行下面的语句产生近似系数cA1、细节系数cD1cA1,cD1 = dwt(s,db1);3. 从系数中构建近似和细节从系数cA1和cD1中构建一层近似A1和细节D1，执行以下代码A1 = upcoef(a,cA1,db1,1,l_s);D1 = upcoef(d,cD1,db1,1,l_s);或A1 = idwt(cA1,db1,l_s);D1 = idwt(,cD1,db1,l_s);4. 显示近似和细节subplot(1,2,1); plot(A1); title(Approximation A1)subplot(1,2,2); plot(D1); title(Detail D1)5. 使用逆小波变换恢复信号A0 = idwt(cA1,cD1,db1,l_s);err = max(abs(s-A0)err =2.2737e-0136. 执行多层小波分解执行3层信号分解C,L = wavedec(s,3,db1);函数返回3层分解的各组分系数C（连接在一个向量里），向量L里返回的是各组分的长度。分解的结构如下7. 抽取近似系数和细节系数从C中抽取3层近似系数cA3 = appcoef(C,L,db1,3);从C中抽取3、2、1层细节系数cD3 = detcoef(C,L,3);cD2 = detcoef(C,L,2);cD1 = detcoef(C,L,1);或者cD1,cD2,cD3 = detcoef(C,L,1,2,3);结果显示如下，从上到下依次为原始信号、3层近似系数和31层细节系数1. 重建3层近似和1、2、3层细节从C中重建3层近似A3 = wrcoef(a,C,L,db1,3);从C中重建1、2、3层细节D1 = wrcoef(d,C,L,db1,1);D2 = wrcoef(d,C,L,db1,2);D3 = wrcoef(d,C,L,db1,3);2. 显示多层分解的结果显示3层分解的结果subplot(2,2,1); plot(A3);title(Approximation A3)subplot(2,2,2); plot(D1);title(Detail D1)subplot(2,2,3); plot(D2);title(Detail D2)subplot(2,2,4); plot(D3);title(Detail D3)10. 从3层分解中重建原始信号A0 = waverec(C,L,db1);err = max(abs(s-A0)err =4.5475e-01311. 粗糙的去噪信号使用小波从信号中移除噪声需要辨识哪个或哪些组分包含噪声，然后重建没有这些组分的信号。在这个例子中，我们注意到连续的近似随着越来越多的高频信息从信号中滤除，噪声变得越来越少。3层近似与原始信号对比会发现变得很干净。对比近似和原始信号，如下subplot(2,1,1);plot(s);title(Original); axis offsubplot(2,1,2);plot(A3);title(Level 3 Approximation);axis off当然，摒弃所有高频信息，我们会失去原始信号中的很多最尖锐的特征。最佳的去噪需要通过一种更精细的叫阈值方法，它只丢弃部分超过一定范围的细节。12. 通过阈值去除噪声先来看3层分析的细节subplot(3,1,1); plot(D1); title(Detail Level 1); axis offsubplot(3,1,2); plot(D2); title(Detail Level 2); axis offsubplot(3,1,3); plot(D3); title(Detail Level 3); axis off从图中可以看到，大多数噪声发生在信号的后面部分，表现在细节上就是出现大波动的地方。如果我们通过设定最大值来限定细节强度，会怎么样呢？这会有降低噪声效果，同时保留不影响的必要细节。这是一种很好的方法。注意到cD1,cD2,cD3是向量，那么我们就可以通过直接操纵这些向量来达到目的，即设置这些向量小于其峰值或平均值的一部分，然后就可以由这些设定了阈值的系数重建新的细节信号D1、D2、D3。实际去噪过程中，可以使用ddencmp函数来计算默认的阈值参数，然后用wdencmp函数来执行实际的去噪过程，代码如下thr,sorh,keepapp = ddencmp(den,wv,s);clean = wdencmp(gbl,C,L,db1,3,thr,sorh,keepapp);注意wdencmp使用了第6步中小波分解的结果C、L，另外指定了db1小波来做分析，指定全局阈值选项gb1.详细参见ddencmp函数和wdencmp函数。显示原始信号和去噪信号如下subplot(2,1,1); plot(s(2000:3920); title(Original)subplot(2,1,2); plot(clean(2000:3920); title(denoised)这里我们只绘制了原始信号中包含噪声的部分，特别注意我们是如何在移除了噪声的情况下仍保持原有的尖锐细节的，这也是小波分析强大的地方。使用命令行函数去除噪声是很笨拙的，而图形工具提供了一种更方便使用的自动化阈值去噪方式。关于去噪过程的更多信息可以在下面的部分找到：* Remove noise from a signal 从一个信号中去除噪声*De-Noising in the Wavelet Toolbox Users Guide 小波工具箱中的除噪*One-Dimensional Variance Adaptive Thresholding of Wavelet Coefficients 一维小波系数自适应阈值*One-Dimensional Variance Adaptive Thresholding of Wavelet Coeffiients in the Wavelet Toolbox Users Guide使用图形接口做一维分析1. 开启一维小波分析工具Wavemenu-Wavelet 1-D2. 加载信号3. 执行一层小波分解使用db1小波执行一层分解4. 放大有关细节5. 执行多层小波分解使用db1小波执行3层分解。选择不同的显示方式：在Display mode下拉菜单下可以选择不同的显示方式，默认的显示方式为Full Decomposition Mode，其他的显示方式及其意义如下Separate Mode：在不同的列中显示细节和近似；Superimpose Mode：在一张图上以不同的颜色显示细节、近似；Tree Mode：显示分解树、原始信号和选择的成分，在分解树上选择你想显示的成分；Show and Scroll Mode：显示3个窗口，第一个显示原始信号和选择的近似信号，第二个显示选择的细节，第三个显示小波系数；Show and Scroll Mode（Stem Cfs）：和Show and Scroll Mode很接近，除了第三个窗口中以杆状图替代颜色条显示小波系数。对于每个分析任务，可以改变默认的显示方式，只要在View-Default Display Mode子菜单下选择理想的方式即可；不同的显示方式会有额外的显示选项，在More Display Options中做选择，这些选项可以控制不同成分的显示、选择是否显示原始信号与细节、近似对比。6. 从信号中移除噪声图形接口提供了以预定义的阈值策略除噪的选项，这使得从信号中除噪非常容易，只需点击De-noise按钮就可以弹出除噪工具。点击Close可以关闭除噪窗口，由于不能同时打开除噪和压缩窗口，所以需要关闭除噪窗口再进行信号压缩。关闭时会提示Update Synthesized Signal提示对话框，点击No，如果点击Yes，合成的信号会加载到主窗口。7. 改善分析图像工具可以在任何时候轻易的改善分析，只需要改变分析的方法就可以了，如使用db3做5层小波分析。8. 压缩信号图形接口提供了自动化或人工阈值的做压缩的功能。默认使用的是全阈值方法，当然也可以使用人工阈值的方法，选择By Level thresholding选项即可，下面的滑动条提供了各级阈值独立调整的功能，相应的调整可以在左边的窗口中看到，在图形窗口中也可以直接拖动来改变阈值。完成选择后，点击压缩按钮可以完成压缩。从压缩的结果可以看到，压缩过程去除了大多数噪声，但保存了信号99.74%的能量。自动化阈值是非常有效的，它使除3.2%的小波系数都归零化了。9. 显示残差点击Residuals按钮可以查看压缩的残差。显示的统计数据包括测量的趋势（平均值、众数、中值）和散布情况（极差、标准差）。另外，工具还提供了概率分布直方图和累计直方图以及时间序列图，如自相关函数、频谱。这些都是和去噪工具是一样的。10. 显示统计分布可以显示一系列有关信号及其组分的统计数据。点击Statistics按钮可以查看统计数据信息，点击Histograms可以查看直方图。从图形接口中导入导出信息保存信息l 保存合成的信号如加载如下信号File Example Analysis Basic Signals with db3 at level 5 Sum of sines，做除噪或压缩处理后，保存合成信号File Save Synthesized Signal，保存后加载文件，会得到如下变量：如果使用除全阈值外的方法时，得到的变量结构如下Synthsig是合成的信号，除噪或压缩的小波方法保存在wname中，相互依赖的各级阈值保存在thrParams中，小波分解的等级数和cell的长度相等，thrParamsi，i从1到5分别保存了阈值间距上下限的值和阈值（间距阈值是允许的，在自适应阈值方法中会用到，参见One-Dimensional Variance Adaptive Thresholding of Wavelet Coefficients）如果使用的全阈值方法，保存的数据结构如下alTHR保存的是全阈值的值。l 保存离散小波变换的系数一个例子的文件内容如下Coefs包含了离散小波变换的系数，longs包含了各组分的长度，thrParams为空，因为合成信号不存在，wname是小波的名字。l 保存分解结果（即保存小波分析的全体数据）小波工具将保存为.wal文件，加载方式为load wdecex1d.wa1 mat文件内容为加载信息加载的文件只要和保存的相应文件中的变量一样即可。Matlab小波工具箱的使用3 (2011-11-11 09:35:01)转载标签：平稳小波分解去噪分类：学科知识一维离散平稳小波分析离散平稳小波分析所用到的函数有swt小波分解和iswt小波重构。在这一部分可以学到如下内容：l 加载信号l 执行平稳小波分解l 从小波系数中构造近似和细节l 显示第一层的近似和细节l 平稳小波逆变换恢复原始信号l 执行多层平稳小波分解l 重构第3层近似l 重构1、2、3层细节l 重构1、2层近似l 显示分解的结果l 从3层分解中重构原始信号l 除噪使用命令行实现一维分析1. 加载信号（一个噪声污染的多普勒效应信号）load noisdopp2. 设置变量s = noisdopp;对于SWT变换，如果需要在第k层分解信号，那么原始信号需要能够平分成2k份。所以如果原始信号的长度不满足要求，需要使用Signal Extension GUI工具或使用wextend函数来扩展它。3. 执行一层平稳小波分解swa,swd = swt(s,1,db1);函数执行将产生1层近似和细节的系数，两者和信号的长度是相等的，这也是平稳小波和普通小波不同的地方，从而使它在某些领域有好的效果。4. 显示近似、细节小波系数显示一层近似和细节的系数subplot(1,2,1), plot(swa); title(Approximation cfs)subplot(1,2,2), plot(swd); title(Detail cfs)5. 使用平稳小波逆变换恢复原始信号A0 = iswt(swa,swd,db1);重构的误差为err = norm(s-A0)err =2.1450e-0146. 从系数中构建近似和细节构建一层近似和细节，输入如下代码nulcfs = zeros(size(swa);A1 = iswt(swa,nulcfs,db1);D1 = iswt(nulcfs,swd,db1);显示结果如下subplot(1,2,1), plot(A1); title(Approximation A1);subplot(1,2,2), plot(D1); title(Detail D1);7. 执行多层平稳小波分解使用db1小波做3层小波分解swa,swd = swt(s,3,db1);产生的近似系数在swa中，细节系数在swd中，而且它们有相同的长度。8. 显示近似和细节的系数kp = 0;for i = 1:3subplot(3,2,kp+1), plot(swa(i,:);title(Approx. cfs level ,num2str(i)subplot(3,2,kp+2), plot(swd(i,:);title(Detail cfs level ,num2str(i)kp = kp + 2;end9. 从系数中重建第3层的近似mzero = zeros(size(swd);A = mzero;A(3,:) = iswt(swa,mzero,db1);10. 从系数中重建细节D = mzero;for i = 1:3swcfs = mzero;swcfs(i,:) = swd(i,:);D(i,:) = iswt(mzero,swcfs,db1);End11. 从第3层近似和第2、3层细节中重建第1、2层近似重建2、3层的近似A(2,:) = A(3,:) + D(3,:);A(1,:) = A(2,:) + D(2,:);显示第1、2、3层的近似和细节kp = 0;for i = 1:3subplot(3,2,kp+1), plot(A(i,:);title(Approx. level ,num2str(i)subplot(3,2,kp+2), plot(D(i,:);title(Detail level ,num2str(i)kp = kp + 2;end12. 阈值除噪要去除噪声，先使用ddencmp函数计算一个默认的全局阈值，再使用wthresh函数执行细节系数的实际阈值，然后使用iswt得到除噪信号。所有在一维离散小波变换中选择阈值的方法在一维平稳小波变换中都是有效的，GUI中两者的使用方法也一样。thr,sorh = ddencmp(den,wv,s);dswd = wthresh(swd,sorh,thr);clean = iswt(swa,dswd,db1);显示原始信号和除噪信号如下subplot(2,1,1), plot(s);title(Original signal)subplot(2,1,2), plot(clean);title(denoised signal)得到的信号仍然有一些噪声，结果可以通过考虑对信号进行5层分解，然后再除噪来改善，如下swa,swd = swt(s,5,db1);thr,sorh = ddencmp(den,wv,s);dswd = wthresh(swd,sorh,thr);clean = iswt(swa,dswd,db1);subplot(2,1,1), plot(s); title(Original signal)subplot(2,1,2), plot(clean); title(denoised signal)第二种实现相同效果的语法如下lev = 5; swc = swt(s,lev,db1);swcden = swc;swcden(1:end-1,:) = wthresh(swcden(1:end-1,:),sorh,thr);clean = iswt(swcden,db1);使用图形接口做一维除噪分析1. 开启一维平稳小波除噪工具输入wavemenu，选择SWT Denoising 1-D，出现如下GUI2. 加载信号Load Signal3. 执行平稳小波分解使用db1小波执行5层小波分解，得到的是非抽取系数（nondecimated coefficients），它们是使用相同的离散小波变换来得到的，只是省略了抽取的步骤。得到的结果如下4. 使用平稳小波变换除噪可以使用GUI默认的参数做除噪处理。右边的滑动条可以控制各级系数的阈值大小，也可以直接在系数图中直接拖动来调整阈值的大小，注意近似系数中没有阈值。点击denoise进行除噪处理得到的效果是非常好的，但似乎在信号不连续的地方出现了过平滑，这个可以从残差图中看出来，在800的位置出现了突降点。选择hard阈值模式代替soft模式，再进行除噪，结果如下可以看到这次效果非常好，而且残差图也看起来像白噪声序列。为了验证这一点，可以点击Residuals按钮查看残差图及相关统计数据来详细说明。离散傅里叶变换DFT基本原理图解 (2011-11-12 11:10:49)转载标签：矩形窗频谱余弦窗函数频域分类：学科知识前两天看书，看到解释DFT基本原理的，认为讲的挺好的，虽然在信号与系统里也学过类似的图，但没有对比，有些东西领会的不深，这个通过对比不同时间窗、不同采样率对频谱泄露的影响等，讲的很好，然后我进一步发挥，对比了矩形窗和汉宁窗，以加深理解。先来看加矩形窗的整周期采样过程（整周期采样体现在窗宽为周期的整数倍），参数如下T0 = 1; % 余弦波的周期f0 = 1/T0; % 余弦波的频率A = 1; % 余弦波的振幅Ts = 0.1; % 时域采样周期Fs = 1/Ts;wT = 1; % 窗的宽度windowflag = 1; % 1为矩形窗，其他为汉宁窗结果如下可以看到整周期采样没有频谱泄露。其实在加窗时就会造成频谱的泄露，如sinc函数形状有明显的旁瓣，但如果窗函数的宽度是信号周期的整数倍，则频域采样（时域周期延拓）将正好采在sinc函数的峰值处，表现出来就是没有频谱泄露（上图中最后一幅图）。如果不是整周期采样，如设窗宽为wT = 1.4; % 窗的宽度结果为可以看到非整周期采样时，频域的采样点采在sinc函数的非峰值处，造成明显的频谱泄露（上图中最后一幅图）。而如果这时使用其他窗函数，则可以改善频谱泄露的状况，如使用汉宁窗windowflag = 2; % 1为矩形窗，其他为汉宁窗结果如下可以看到汉宁窗的旁瓣非常低，能量相对矩形窗非常集中，最后的结果是虽然仍有频谱泄露，但相对于矩形窗已经好得多了，主要还在于其压低的旁瓣很小。可见加窗的效果。汉宁窗也有它的缺点，就是它的峰值处的带宽比矩形窗要宽，可以说分辨率没有矩形窗好。注意到汉宁窗频谱的高度只有不到1，而矩形窗的频谱则达1.5，这个没有关系，因为汉宁窗在时域压低了边缘的幅度，相当于减小了信号的能量，所以频域的峰值会降低，这是可以也需要在频域处理完反变换回来时除以汉宁窗来解除其压低能量的效果的。记得学数字信号处理时，说造成频谱泄露的两个原因是窗函数有能量泄露（上面如实展示了）、周期延拓造成不连续点。对于后者，个人认为值得商榷，造成不连续点只是表象，如果说这是本质原因，那采样值总是离散值，离散值就是不连续点，总会造成频谱泄露的呀，但对整周期采样的余弦波，虽然是离散点，但没有频谱泄露（可能有，但很小），所以其真正原因是没有整周期采样，但问题又来了，实际中的信号可不是简单的余弦波、正弦波，有固定频率的，所以无法整周期采样，整谁的周期啊？所以频谱泄露只能抑制，不能消除。以上结果的代码如下main.m文件clcclear allclose allts = 0.001; % 连续函数的显示采样周期fs = 0.01; % 频域内的显示采样周期T0 = 1; % 余弦波的周期f0 = 1/T0; % 余弦波的频率A = 1; % 余弦波的振幅Ts = 0.1; % 时域采样周期Fs = 1/Ts;wT = 1.4; % 窗的宽度windowflag = 2; % 1为矩形窗，其他为汉宁窗handle.txrange = -T0 T0*2; % 时域的显示范围handle.tyrange = -A*1.5 A*1.5;handle.fxrange = -f0 f0*13; % 频域的显示范围handle.fyrange = -A*0.5 A*2;% 余弦信号t = handle.txrange(1):ts:handle.txrange(2);xt = A*cos(2*pi*f0*t);index = 1;handle.data(index).x = t;handle.data(index).y = xt;handle.data(index).type = plot;handle.data(index).ylim = -A A*1.2;handle.data(index).title = 余弦信号x(t);% 余弦信号的频谱f = handle.fxrange(1):fs:handle.fxrange(2);xf = A*(f = f0 | f = -f0);index = 2;handle.data(index).x = f;handle.data(index).y = xf;handle.data(index).type = stem;handle.data(index).ylim = -A*0.2 A*1.2;handle.data(index).title = 余弦信号频谱x(f);% 采样序列（周期脉冲序列）t = handle.txrange(1):Ts:handle.txrange(2);pt = ones(size(t);index = 3;handle.data(index).x = t;handle.data(index).y = pt;handle.data(index).type = stem;handle.data(index).ylim = -A*0.2 A*1.2;handle.data(index).title = 采样脉冲p(t);% 采样序列的频谱sampletemp = round(handle.fxrange(1)/fs):round(handle.fxrange(2)/fs);f = sampletemp*fs;pf = mod(sampletemp, round(Fs/fs) = 0;index = 4;handle.data(index).x = f;handle.data(index).y = pf;handle.data(index).type = stem;handle.data(index).ylim = -A*0.2 A*1.2;handle.data(index).title = 采样脉冲频谱p(f);% 采样余弦波t = handle.txrange(1):Ts:handle.txrange(2);xpt = A*cos(2*pi*f0*t);index = 5;handle.data(index).x = t;handle.data(index).y = xpt;handle.data(index).type = stem;handle.data(index).ylim = -A A*1.2;handle.data(index).title = 采样信号x(t)p(t);% 采样余弦波的频谱f = handle.fxrange(1):fs:handle.fxrange(2);xpf = conv(xf, double(pf), same);index = 6;handle.data(index).x = f;handle.data(index).y = xpf;handle.data(index).type = stem;handle.data(index).ylim = -A*0.2 A*1.2;handle.data(index).title = 采样信号频谱x(f)*p(f);% 窗函数t = handle.txrange(1):ts:handle.txrange(2);if windowflag = 1wt = t=0 & t=wT;elsewt = sin(t*pi/wT).2;wt(twT) = 0;endindex = 7;handle.data(index).x = t;handle.data(index).y = wt;handle.data(index).type = plot;handle.data(index).ylim = -A*0.2 A*1.2;handle.data(index).title = 窗函数w(t);% 窗函数的频谱f = handle.fxrange(1):fs:handle.fxrange(2);if windowflag = 1wf = sinc(wT*f)*wT; % 矩形窗elsewf = wT/2*sinc(wT*f) + wT/4*sinc(wT*f-1) + wT/4*sinc(wT*f-1); % 汉宁窗endindex = 8;handle.data(index).x = f;handle.data(index).y = abs(wf);handle.data(index).type = plot;handle.data(index).ylim = -A*0.2 A*wT*1.2;handle.data(index).title = 窗函数频谱|w(f)|;% 采样余弦信号加窗t = 0:Ts:wT-Ts/100;if windowflag = 1xpwt = A.*cos(2*pi*f0*t);elsexpwt = sin(t*pi/wT).2*A.*cos(2*pi*f0*t);endindex = 9;handle.data(index).x = t;handle.data(index).y = xpwt;handle.data(index).type = stem;handle.data(index).ylim = -A A*1.2;handle.data(index).title = 信号加窗x(t)p(t)w(t);% 加窗后的频谱f = handle.fxrange(1):fs:handle.fxrange(2);xpwf = conv(xpf, wf, same);index = 10;handle.data(index).x = f;handle.data(index).y = abs(xpwf);handle.data(index).type = plot;handle.data(index).ylim = -A*0.2 A*wT*1.2;handle.data(index).title = 加窗后频谱|x(f)*p(f)*w(f)|