Matlab频谱分析程序
利用Matlab绘制正弦信号的频谱图并做相关分析
利用Matlab绘制正弦信号的频谱图并做相关分析一、作业要求：1、信号可变（信号的赋值、相位、频率可变）；2、采样频率fs可变；3、加各种不同的窗函数并分析其影响；4、频谱校正；5、频谱细化。
二、采用matlab编写如下程序：clear;clf;fs=100;N=1024; %采样频率和数据点数A=20;B=30;C=0.38;n=0:N-1;t=n/fs; %时间序列x=A*sin(2*pi*B*t+C); %信号y=fft(x,N); %对信号进行傅里叶变换yy=abs(y); %求得傅里叶变换后的振幅yy=yy*2/N; %幅值处理f=n*fs/N; %频率序列subplot(3,3,1),plot(f,yy); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图1：fs=100，N=1024');grid on;%两种信号叠加，x=A*sin(2*pi*B*t+C)+2*A*sin(2*pi*1.5*B*t+2.5*C); %信号y=fft(x,N); %对信号进行傅里叶变换yy=abs(y); %求得傅里叶变换后的振幅yy=yy*2/N; %幅值处理f=n*fs/N; %频率序列subplot(3,3,2),plot(f,yy); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图2：fs=100,N=1024，两种信号叠加');grid on;%加噪声之后的图像x=A*sin(2*pi*B*t+C)+28*randn(size(t));y=fft(x,N);yy=abs(y);yy=yy*2/N; %幅值处理subplot(3,3,3),plot(f(1:N/2.56),yy(1:N/2.56));xlabel('频率/\itHz');ylabel('振幅');title('图3：fs=100,N=1024混入噪声');grid on;%改变采样点数N=128N=128;n=0:N-1;t=n/fs; %时间序列x=A*sin(2*pi*B*t+C); %信号y=fft(x,N); %对信号进行傅里叶变换yy=abs(y); %求得傅里叶变换后的振幅yy=yy*2/N; %幅值处理f=n*fs/N; %频率序列subplot(3,3,4),plot(f(1:N/2.56),yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图4：fs=100，N=128');grid on;%改变采样频率为200Hz时的频谱fs=400;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信号y=fft(x,N); %对信号进行快速傅里叶变换yy=abs(y); %求取傅里叶变换的振幅yy=yy*2/N; %幅值处理f=n*fs/N;subplot(3,3,5),plot(f(1:N/2.56),yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图5：fs=400，N=1024');grid on;%加三角窗函数fs=100;N=1024; %采样频率和数据点数n=0:N-1;t=n/fs; %时间序列x=A*sin(2*pi*B*t+C); %信号window=triang(N);%生成三角窗函数x=x.*window';%加窗函数y=fft(x,N); %对信号进行傅里叶变换yy=abs(y); %求得傅里叶变换后的振幅yy=yy*2/N; %幅值处理f=n*fs/N; %频率序列subplot(3,3,6),plot(f(1:N/2.56),2*yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图6：fs=100,N=1024,加三角窗函数');grid on;%加海明窗函数后的频谱fs=100;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信号window=hamming(N);%生成海明窗函数x=x.*window';%加窗函数y=fft(x,N); %对信号进行快速傅里叶变换yy=abs(y); %求取傅里叶变换的振幅yy=yy*2/N; %幅值处理f=n*fs/N;subplot(3,3,7),plot(f(1:N/2.56),1.852*yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图7：fs=100，N=1024,加海明窗函数');grid on;%加汉宁窗函数后的频谱fs=100;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信号window=hanning(N);%生成汉宁窗函数x=x.*window';%加窗函数y=fft(x,N); %对信号进行快速傅里叶变换yy=abs(y); %求取傅里叶变换的振幅yy=yy*2/N; %幅值处理f=n*fs/N;subplot(3,3,8),plot(f(1:N/2.56),2*yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图8：fs=100，N=1024,加汉宁窗函数');grid on;三、运行结果如下：四、分析与结论：1）从所做图像可以看出，信号的幅值均小于真实值，说明在截断信号时存在泄露。
实验2利用MATLAB分析信号频谱及系统的频率特性
实验2利用MATLAB分析信号频谱及系统的频率特性引言：在信号处理和通信领域中，频谱分析是一项非常重要的技术。
频谱分析可以帮助我们了解信号的频率特性，包括频率成分和幅度。
MATLAB是一款功能强大的数学软件，提供了多种工具和函数用于信号处理和频谱分析。
本实验旨在通过MATLAB分析信号频谱及系统的频率特性，深入理解信号处理和频域分析的原理和应用。
实验步骤：1.生成一个信号并绘制其时域波形。
首先，我们可以使用MATLAB提供的函数生成一个信号。
例如，我们可以生成一个用正弦函数表示的周期信号。
```matlabt=0:0.001:1;%时间范围为0到1秒，采样率为1000Hzf=10;%信号频率为10Hzx = sin(2*pi*f*t); % 生成正弦信号plot(t,x) % 绘制信号的时域波形图title('Time domain waveform') % 添加标题```2.计算信号的频谱并绘制频谱图。
使用MATLAB中的FFT函数可以计算信号的频谱。
FFT函数将信号从时域转换为频域。
```matlabFs=1000;%采样率为1000HzL = length(x); % 信号长度NFFT = 2^nextpow2(L); % FFT长度X = fft(x,NFFT)/L; % 计算X(k)f = Fs/2*linspace(0,1,NFFT/2+1); % 计算频率轴plot(f,2*abs(X(1:NFFT/2+1))) % 绘制频谱图title('Frequency spectrum') % 添加标题```3.使用MATLAB分析系统的频率特性。
MATLAB提供了Signal Processing Toolbox，其中包含了分析系统频率特性的函数和工具。
```matlabHd = designfilt('lowpassfir', 'FilterOrder', 6,'CutoffFrequency', 0.3, 'SampleRate', Fs); % 设计一个低通滤波器fvtool(Hd) % 显示滤波器的频率响应``````matlab[W,F] = freqz(Hd); % 计算滤波器的频率响应plot(F,abs(W)) % 绘制滤波器的振幅响应title('Frequency response of lowpass filter') % 添加标题```实验结果：运行上述代码后，我们可以得到如下结果：1.时域波形图2.频谱图3.滤波器频率响应讨论与结论：本实验通过MATLAB分析信号频谱及系统的频率特性，深入理解了信号处理和频域分析的原理和应用。
利用Matlab进行频谱分析的方法
利用Matlab进行频谱分析的方法引言频谱分析是信号处理和电子工程领域中一项重要的技术，用于分析信号在频率域上的特征和频率成分。
在实际应用中，频谱分析广泛应用于音频处理、图像处理、通信系统等领域。
Matlab是一种强大的工具，可以提供许多功能用于频谱分析。
本文将介绍利用Matlab进行频谱分析的方法和一些常用的工具。
一、Matlab中的FFT函数Matlab中的FFT（快速傅里叶变换）函数是一种常用的频谱分析工具。
通过使用FFT函数，我们可以将时域信号转换为频域信号，并得到信号的频谱特征。
FFT 函数的使用方法如下：```Y = fft(X);```其中，X是输入信号，Y是输出的频域信号。
通过该函数，我们可以得到输入信号的幅度谱和相位谱。
二、频谱图的绘制在进行频谱分析时，频谱图是一种直观和易于理解的展示形式。
Matlab中可以使用plot函数绘制频谱图。
首先，我们需要获取频域信号的幅度谱。
然后，使用plot函数将频率与幅度谱进行绘制。
下面是一个示例：```X = 1:1000; % 时间序列Y = sin(2*pi*10*X) + sin(2*pi*50*X); % 输入信号Fs = 1000; % 采样率N = length(Y); % 信号长度Y_FFT = abs(fft(Y)); % 计算频域信号的幅度谱f = (0:N-1)*(Fs/N); % 频率坐标plot(f, Y_FFT);```通过上述代码，我们可以得到输入信号在频谱上的特征，并将其可视化为频谱图。
三、频谱分析的应用举例频谱分析可以应用于许多实际问题中。
下面将介绍两个常见的应用举例：语音信号分析和图像处理。
1. 语音信号分析语音信号分析是频谱分析的一个重要应用领域。
通过对语音信号进行频谱分析，我们可以探索声波的频率特性和信号的频率成分。
在Matlab中，可以使用wavread 函数读取音频文件，并进行频谱分析。
下面是一个示例：```[waveform, Fs] = wavread('speech.wav'); % 读取音频文件N = length(waveform); % 信号长度waveform_FFT = abs(fft(waveform)); % 计算频域信号的幅度谱f = (0:N-1)*(Fs/N); % 频率坐标plot(f, waveform_FFT);```通过上述代码，我们可以获取语音信号的频谱特征，并将其可视化为频谱图。
利用Matlab绘制正弦信号的频谱图并做相关分析
利用Matlab绘制正弦信号的频谱图并做相关分析一、作业要求：1、信号可变（信号的赋值、相位、频率可变）；2、采样频率fs可变；3、加各种不同的窗函数并分析其影响；4、频谱校正；5、频谱细化。
二、采用matlab编写如下程序：clear;clf;fs=100;N=1024; %采样频率和数据点数A=20;B=30;C=0.38;n=0:N-1;t=n/fs; %时间序列x=A*sin(2*pi*B*t+C); %信号y=fft(x,N); %对信号进行傅里叶变换yy=abs(y); %求得傅里叶变换后的振幅yy=yy*2/N; %幅值处理f=n*fs/N; %频率序列subplot(3,3,1),plot(f,yy); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图1：fs=100，N=1024');grid on;%两种信号叠加，x=A*sin(2*pi*B*t+C)+2*A*sin(2*pi*1.5*B*t+2.5*C); %信号y=fft(x,N); %对信号进行傅里叶变换yy=abs(y); %求得傅里叶变换后的振幅yy=yy*2/N; %幅值处理f=n*fs/N; %频率序列subplot(3,3,2),plot(f,yy); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图2：fs=100,N=1024，两种信号叠加');grid on;%加噪声之后的图像x=A*sin(2*pi*B*t+C)+28*randn(size(t));y=fft(x,N);yy=abs(y);yy=yy*2/N; %幅值处理subplot(3,3,3),plot(f(1:N/2.56),yy(1:N/2.56));xlabel('频率/\itHz');ylabel('振幅');title('图3：fs=100,N=1024混入噪声');grid on;%改变采样点数N=128N=128;n=0:N-1;t=n/fs; %时间序列x=A*sin(2*pi*B*t+C); %信号y=fft(x,N); %对信号进行傅里叶变换yy=abs(y); %求得傅里叶变换后的振幅yy=yy*2/N; %幅值处理f=n*fs/N; %频率序列subplot(3,3,4),plot(f(1:N/2.56),yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图4：fs=100，N=128');grid on;%改变采样频率为200Hz时的频谱fs=400;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信号y=fft(x,N); %对信号进行快速傅里叶变换yy=abs(y); %求取傅里叶变换的振幅yy=yy*2/N; %幅值处理f=n*fs/N;subplot(3,3,5),plot(f(1:N/2.56),yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图5：fs=400，N=1024');grid on;%加三角窗函数fs=100;N=1024; %采样频率和数据点数n=0:N-1;t=n/fs; %时间序列x=A*sin(2*pi*B*t+C); %信号window=triang(N);%生成三角窗函数x=x.*window';%加窗函数y=fft(x,N); %对信号进行傅里叶变换yy=abs(y); %求得傅里叶变换后的振幅yy=yy*2/N; %幅值处理f=n*fs/N; %频率序列subplot(3,3,6),plot(f(1:N/2.56),2*yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图6：fs=100,N=1024,加三角窗函数');grid on;%加海明窗函数后的频谱fs=100;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信号window=hamming(N);%生成海明窗函数x=x.*window';%加窗函数y=fft(x,N); %对信号进行快速傅里叶变换yy=abs(y); %求取傅里叶变换的振幅yy=yy*2/N; %幅值处理f=n*fs/N;subplot(3,3,7),plot(f(1:N/2.56),1.852*yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图7：fs=100，N=1024,加海明窗函数');grid on;%加汉宁窗函数后的频谱fs=100;N=1024;n=0:N-1;t=n/fs;x=A*sin(2*pi*B*t+C); %信号window=hanning(N);%生成汉宁窗函数x=x.*window';%加窗函数y=fft(x,N); %对信号进行快速傅里叶变换yy=abs(y); %求取傅里叶变换的振幅yy=yy*2/N; %幅值处理f=n*fs/N;subplot(3,3,8),plot(f(1:N/2.56),2*yy(1:N/2.56)); %绘出随频率变化的振幅xlabel('频率/\itHz');ylabel('振幅');title('图8：fs=100，N=1024,加汉宁窗函数');grid on;三、运行结果如下：四、分析与结论：1）从所做图像可以看出，信号的幅值均小于真实值，说明在截断信号时存在泄露。
实验三用FFT对信号进行频谱分析和MATLAB程序
实验三用FFT对信号进行频谱分析和MATLAB程序实验三中使用FFT对信号进行频谱分析的目的是通过将时域信号转换为频域信号，来获取信号的频谱信息。
MATLAB提供了方便易用的函数来实现FFT。
首先，我们需要了解FFT的原理。
FFT（快速傅里叶变换）是一种快速计算离散傅里叶变换（DFT）的算法，用于将离散的时间域信号转换为连续的频域信号。
FFT算法的主要思想是将问题划分为多个规模较小的子问题，并利用DFT的对称性质进行递归计算。
FFT算法能够帮助我们高效地进行频谱分析。
下面是一个使用MATLAB进行频谱分析的示例程序：```matlab%生成一个10秒钟的正弦波信号，频率为1Hz，采样率为100Hzfs = 100; % 采样率t = 0:1/fs:10-1/fs; % 时间范围f=1;%正弦波频率x = sin(2*pi*f*t);%进行FFT计算N = length(x); % 信号长度X = fft(x); % FFT计算magX = abs(X)/N; % 幅值谱frequencies = (0:N-1)*(fs/N); % 频率范围%绘制频谱图figure;plot(frequencies, magX);xlabel('频率(Hz)');ylabel('振幅');title('信号频谱');```上述代码生成了一个10秒钟的正弦波信号，频率为1 Hz，采样率为100 Hz。
通过调用MATLAB的fft函数计算信号的FFT，然后计算每个频率分量的幅值谱，并绘制出信号频谱图。
在频谱图中，横轴表示频率，纵轴表示振幅。
该实验需要注意以下几点：1.信号的采样率要与信号中最高频率成一定比例，以避免采样率不足导致的伪频谱。
2.FFT计算结果是一个复数数组，我们一般只关注其幅值谱。
3.频率范围是0到采样率之间的频率。
实验三的报告可以包含以下内容：1.实验目的和背景介绍。
基于Matlab的DFT及FFT频谱分析
基于Matlab的DFT及FFT频谱分析基于Matlab的DFT及FFT频谱分析一、引言频谱分析是信号处理中的重要任务之一，它可以揭示信号的频率特性和能量分布。
离散傅里叶变换（DFT）及快速傅里叶变换（FFT）是常用的频谱分析工具，广泛应用于许多领域。
本文将介绍通过Matlab进行DFT及FFT频谱分析的方法和步骤，并以实例详细说明。
二、DFT及FFT原理DFT是一种将时域信号转换为频域信号的离散变换方法。
它将信号分解成若干个正弦和余弦函数的叠加，得到频率和幅度信息。
FFT是一种高效的计算DFT的算法，它利用信号的对称性和周期性，将计算复杂度从O(N^2)降低到O(NlogN)。
FFT通过将信号分解成不同长度的子序列，递归地进行计算，最终得到频谱信息。
三、Matlab中的DFT及FFT函数在Matlab中，DFT及FFT可以通过内置函数进行计算。
其中，DFT使用函数fft，FFT使用函数fftshift。
fft函数可直接计算信号的频谱，fftshift函数对频谱进行频移操作，将低频移到频谱中心。
四、Matlab中DFT及FFT频谱分析步骤1. 读取信号数据首先，将待分析的信号数据读入到Matlab中。
可以使用内置函数load读取文本文件中的数据，或通过自定义函数生成模拟信号数据。
2. 时域分析通过plot函数将信号数据在时域进行绘制，以观察信号的波形。
可以设置合适的坐标轴范围和标签，使图像更加清晰。
3. 信号预处理针对不同的信号特点，可以进行预处理操作，例如去除直流分量、滤波等。
这些操作可提高信号的频谱分析效果。
4. 计算DFT/FFT使用fft函数计算信号数据的DFT/FFT，并得到频谱。
将信号数据作为输入参数，设置采样频率和点数，计算得到频谱数据。
5. 频域分析通过plot函数将频谱数据在频域进行绘制，观察信号的频率特性。
可以设置合适的坐标轴范围和标签，使图像更加清晰。
6. 结果解读根据频谱图像，分析信号的频率成分、幅度分布和峰值位置。
MATLAB中FFT的使用方法(频谱分析)
说明：以下资源来源于《数字信号处理的MATLAB实现》万永革主编一.调用方法X=FFT(x)；X=FFT(x，N)；x=IFFT(X);x=IFFT(X,N)用MATLAB进行谱分析时注意：（1）函数FFT返回值的数据结构具有对称性。
例：N=8;n=0:N-1;xn=[4 3 2 6 7 8 9 0];Xk=fft(xn)→Xk =39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929iXk与xn的维数相同，共有8个元素。
Xk的第一个数对应于直流分量，即频率值为0。
（2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。
在IFFT时已经做了处理。
要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。
二.FFT应用举例例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。
采样频率fs=100Hz，分别绘制N=128、1024点幅频图。
clf;fs=100;N=128; %采样频率和数据点数n=0:N-1;t=n/fs; %时间序列x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N); %对信号进行快速Fourier变换mag=abs(y); %求得Fourier变换后的振幅f=n*fs/N; %频率序列subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=128');grid on;%对信号采样数据为1024点的处理fs=100;N=1024;n=0:N-1;t=n/fs;x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号y=fft(x,N); %对信号进行快速Fourier变换mag=abs(y); %求取Fourier变换的振幅f=n*fs/N;subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=1024');grid on;subplot(2,2,4)plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz');ylabel('振幅');title('N=1024');grid on;运行结果：fs=100Hz，Nyquist频率为fs/2=50Hz。
用MATLAB对信号做频谱分析
⽤MATLAB对信号做频谱分析1.⾸先学习下傅⾥叶变换的东西。
学⾼数的时候⽼师只是将傅⾥叶变换简单的说了下，并没有深⼊的讲解。
⽽现在看来，傅⾥叶变换似乎是信号处理的⽅⾯的重点只是呢，现在就先学习学习傅⾥叶变换吧。
上⾯这幅图在知乎⼀个很著名的关于傅⾥叶变换的⽂章中的核⼼插图，我觉得这幅图很直观的就说明了傅⾥叶变换的实质。
时域上的东西直观的反应到了频域上了，很完美的结合到了⼀起，233333. ⽆数正弦波叠加，震荡的叠加的最后结果竟然是⽅波，同理，任何周期性函数竟然都能拆分为傅⾥叶级数的形式，这样的简介与优雅，真令⼈折服。
2.MATLAB对信号做频谱分析代码：(1)对 f1 = Sa(2t)的频谱分析1 clear;clc;2 hold on;3 R=0.05;4 t=-1.2:R:1.2;5 t1 = 2*t;6 f1=sinc(t1); %Sa函数7 subplot(1,2,1),plot(t,f1)8 xlabel('t'),ylabel('f1')9 axis([-2,2,-0.3,1.2]); %写出Sa函数上下限1011 N=1000;12 k=-N:N;13 W1=40;14 W=k*W1/N;15 F=f1*exp(-j*t'*W)*R; %f1的傅⾥叶变换16 F=real(F); %取F的实部17 subplot(1,2,2),plot(W,F)18 xlabel('W'),ylabel('F(jw)')View Code结果如下图：(2)对 f2 = u(t+2) - u(t-2)的频谱分析1 R=0.05;2 t=-3:R:3;3 f2=(t>=-2)-(t>=2);4 subplot(1,2,1),plot(t,f2)5 grid on;6 xlabel('t'),ylabel('f2')7 axis([-3,3,-0.5,1.5]);89 N=1000;k=-N:N;10 W1=40;11 W=k*W1/N;12 F=f2*exp(-j*t'*W)*R;13 F=real(F);14 subplot(1,2,2),plot(W,F)15 grid on;16 xlabel('W'),ylabel('F(jw)')View Code结果如下图：(3)对f3 = t[u(t+1) - u(t-1) ]的频谱分析1 R=0.05;2 h=0.001;3 t=-1.2:R:1.2;4 y=t.*(t>=-1)-t.*(t>=1);5 f4=diff(y)/h;6 subplot(1,2,1),plot(t,y)7 xlabel('t'),ylabel('y')8 axis([-1.2,1.2,-1.2,1.2]);910 N=1000;11 k=-N:N;12 W1=40;13 W=k*W1/N;14 F=y*exp(-j*t'*W)*R;15 F=real(F);16 subplot(1,2,2),plot(W,F)17 xlabel('W'),ylabel('F(jw)')18 axis([-40,40,-0.06,0.06]);View Code结果如下图：(4)对正弦波做FFT频谱分析1 %*************************************************************************%2 % FFT实践及频谱分析 %3 %*************************************************************************%4 %***************正弦波****************%5 fs=100;%设定采样频率6 N=128;7 n=0:N-1;8 t=n/fs;9 f0=10;%设定正弦信号频率10 %⽣成正弦信号11 x=sin(2*pi*f0*t);12 figure(1);13 subplot(231);14 plot(t,x);%作正弦信号的时域波形15 xlabel('t');16 ylabel('y');17 title('正弦信号y=2*pi*10t时域波形');18 grid;1920 %进⾏FFT变换并做频谱图21 y=fft(x,N);%进⾏fft变换22 mag=abs(y);%求幅值23 f=(0:length(y)-1)'*fs/length(y);%进⾏对应的频率转换24 figure(1);25 subplot(232);26 plot(f,mag);%做频谱图27 axis([0,100,0,80]);28 xlabel('频率(Hz)');29 ylabel('幅值');30 title('正弦信号y=2*pi*10t幅频谱图N=128');31 grid;3233 %求均⽅根谱34 sq=abs(y);35 figure(1);36 subplot(233);37 plot(f,sq);38 xlabel('频率(Hz)');39 ylabel('均⽅根谱');40 title('正弦信号y=2*pi*10t均⽅根谱');41 grid;4243 %求功率谱44 power=sq.^2;45 figure(1);46 subplot(234);47 plot(f,power);48 xlabel('频率(Hz)');49 ylabel('功率谱');50 title('正弦信号y=2*pi*10t功率谱');51 grid;5253 %求对数谱54 ln=log(sq);55 figure(1);56 subplot(235);57 plot(f,ln);58 xlabel('频率(Hz)');59 ylabel('对数谱');60 title('正弦信号y=2*pi*10t对数谱');61 grid;6263 %⽤IFFT恢复原始信号64 xifft=ifft(y);65 magx=real(xifft);66 ti=[0:length(xifft)-1]/fs;67 figure(1);68 subplot(236);69 plot(ti,magx);70 xlabel('t');71 ylabel('y');72 title('通过IFFT转换的正弦信号波形');73 grid;View Code执⾏结果如下图：(5)对矩形波做FFT频谱分析1 %****************2.矩形波****************%2 fs=10;%设定采样频率3 t=-5:0.1:5;4 x=rectpuls(t,2);5 x=x(1:99);6 figure(1);7 subplot(231); plot(t(1:99),x);%作矩形波的时域波形8 xlabel('t');9 ylabel('y');10 title('矩形波时域波形');11 grid;1213 %进⾏FFT变换并做频谱图14 y=fft(x);%进⾏fft变换15 mag=abs(y);%求幅值16 f=(0:length(y)-1)'*fs/length(y);%进⾏对应的频率转换17 figure(1);18 subplot(232);19 plot(f,mag);%做频谱图20 xlabel('频率(Hz)');21 ylabel('幅值');22 title('矩形波幅频谱图');23 grid;2425 %求均⽅根谱26 sq=abs(y);27 figure(1);28 subplot(233);29 plot(f,sq);30 xlabel('频率(Hz)');31 ylabel('均⽅根谱');32 title('矩形波均⽅根谱');33 grid;3435 %求功率谱36 power=sq.^2;37 figure(1);38 subplot(234);39 plot(f,power);40 xlabel('频率(Hz)');41 ylabel('功率谱');42 title('矩形波功率谱');43 grid;4445 %求对数谱46 ln=log(sq);47 figure(1);48 subplot(235);49 plot(f,ln);50 xlabel('频率(Hz)');51 ylabel('对数谱');52 title('矩形波对数谱');53 grid;5455 %⽤IFFT恢复原始信号56 xifft=ifft(y);57 magx=real(xifft);58 ti=[0:length(xifft)-1]/fs;59 figure(1);60 subplot(236);61 plot(ti,magx);62 xlabel('t');63 ylabel('y');64 title('通过IFFT转换的矩形波波形');65 grid;View Code执⾏结果如下图：(6)对⽩噪声做频谱分析1 %****************3.⽩噪声****************%2 fs=10;%设定采样频率3 t=-5:0.1:5;4 x=zeros(1,100);5 x(50)=100000;6 figure(1);7 subplot(231);8 plot(t(1:100),x);%作⽩噪声的时域波形9 xlabel('t');10 ylabel('y');11 title('⽩噪声时域波形');12 grid;1314 %进⾏FFT变换并做频谱图15 y=fft(x); %进⾏fft变换16 mag=abs(y);%求幅值17 f=(0:length(y)-1)'*fs/length(y);%进⾏对应的频率转换18 figure(1);19 subplot(232);20 plot(f,mag);%做频谱图21 xlabel('频率(Hz)');22 ylabel('幅值');23 title('⽩噪声幅频谱图');24 grid;2526 %求均⽅根谱27 sq=abs(y);28 figure(1);29 subplot(233);30 plot(f,sq);31 xlabel('频率(Hz)');32 ylabel('均⽅根谱');33 title('⽩噪声均⽅根谱');34 grid;3536 %求功率谱37 power=sq.^2;38 figure(1);39 subplot(234);40 plot(f,power);41 xlabel('频率(Hz)');42 ylabel('功率谱');43 title('⽩噪声功率谱');44 grid;4546 %求对数谱47 ln=log(sq);48 figure(1);49 subplot(235);50 plot(f,ln);51 xlabel('频率(Hz)');52 ylabel('对数谱');53 title('⽩噪声对数谱');54 grid;5556 %⽤IFFT恢复原始信号57 xifft=ifft(y);58 magx=real(xifft);59 ti=[0:length(xifft)-1]/fs;60 figure(1);61 subplot(236);62 plot(ti,magx);63 xlabel('t');64 ylabel('y');65 title('通过IFFT转换的⽩噪声波形');66 grid;View Code执⾏结果如下：。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Matlab 信号处理工具箱 谱估计专题频谱分析Spectral estimation （谱估计）的目标是基于一个有限的数据集合描述一个信号的功率（在频率上的）分布。
功率谱估计在很多场合下都是有用的，包括对宽带噪声湮没下的信号的检测。
从数学上看，一个平稳随机过程n x 的power spectrum （功率谱）和correlation sequence （相关序列）通过discrete-time Fourier transform （离散时间傅立叶变换）构成联系。
从normalized frequency （归一化角频率）角度看，有下式()()j mxx xxm S R m eωω∞-=-∞=∑注：()()2xx S X ωω=，其中()/2/2limN j n n N n N X x e ωω=-=∑πωπ-<≤。
其matlab近似为X=fft(x,N)/sqrt(N)，在下文中()L X f 就是指matlab fft 函数的计算结果了使用关系2/s f f ωπ=可以写成物理频率f 的函数，其中s f 是采样频率()()2/sjfm f xx xxm S f R m eπ∞-=-∞=∑相关序列可以从功率谱用IDFT 变换求得：()()()/22//22sss f jfm f j m xx xx xx sf S e S f e R m d df f πωππωωπ--==⎰⎰序列n x 在整个Nyquist 间隔上的平均功率可以表示为()()()/2/202ss f xx xx xx sf S S f R d df f ππωωπ--==⎰⎰上式中的()()2xx xx S P ωωπ=以及()()xx xx sS f P f f = 被定义为平稳随机信号n x 的power spectral density (PSD)（功率谱密度） 一个信号在频带[]1212,,0ωωωωπ≤<≤上的平均功率可以通过对PSD 在频带上积分求出[]()()211212,xxxxP P d P d ωωωωωωωωωω--=+⎰⎰ 从上式中可以看出()xx P ω是一个信号在一个无穷小频带上的功率浓度，这也是为什么它叫做功率谱密度。
PSD 的单位是功率（ 瓦特）每单位频率。
在()xx P ω的情况下，这是瓦特/弧度/抽或只是瓦特/弧度。
在()xx P f 的情况下单位是瓦特/赫兹。
PSD 对频率的积分得到的单位是瓦特，正如平均功率[]12,P ωω所期望的那样。
对实信号，PSD 是关于直流信号对称的，所以0ωπ≤≤的()xx P ω就足够完整的描述PSD 了。
然而要获得整个Nyquist 间隔上的平均功率，有必要引入单边PSD 的概念：()()0020onesided xxP P πωωωωπ-≤<⎧=⎨≤<⎩信号在频带[]1212,,0ωωωωπ≤<≤上的平均功率可以用单边PSD 求出[]()2121,onesidedP P d ωωωωωω=⎰频谱估计方法Matlab 信号处理工具箱提供了三种方法PSD直接从信号本身估计出来。
最简单的就是periodogram（周期图法），一种改进的周期图法是Welch's method。
更现代的一种方法是multitaper method（多椎体法）。
Parametric methods （参量类方法）这类方法是假设信号是一个由白噪声驱动的线性系统的输出。
这类方法的例子是Yule-Walker autoregressive (AR) method和Burg method。
这些方法先估计假设的产生信号的线性系统的参数。
这些方法想要对可用数据相对较少的情况产生优于传统非参数方法的结果。
Subspace methods （子空间类）又称为high-resolution methods（高分辨率法）或者super-resolution methods（超分辨率方法）基于对自相关矩阵的特征分析或者特征值分解产生信号的频率分量。
代表方法有multiple signal classification (MUSIC) method或eigenvector (EV) method。
这类方法对线谱（正弦信号的谱）最合适，对检测噪声下的正弦信号很有效，特别是低信噪比的情况。
Nonparametric Methods非参数法下面讨论periodogram, modified periodogram, Welch, 和 multitaper法。
同时也讨论CPSD函数，传输函数估计和相关函数。
Periodogram周期图法一个估计功率谱的简单方法是直接求随机过程抽样的DFT，然后取结果的幅度的平方。
这样的方法叫做周期图法。
x n的PSD的周期图估计是一个长L的信号[]L注：这里()L X f 运用的是matlab 里面的fft 的定义不带归一化系数，所以要除以L 其中()[]12/0s L jfn f L L n X f x n e π--==∑实际对()L X f 的计算可以只在有限的频率点上执行并且使用FFT 。
实践上大多数周期图法的应用都计算N 点PSD 估计0,1,,1N -其中()[]12/0L jkn N L k L n X f x n e π--==∑选择N 是大于L 的下一个2的幂次是明智的，要计算[]L k X f 我们直接对[]L x n 补零到长度为N 。
假如L>N ，在计算[]L k X f 前，我们必须绕回[]L x n 模N 。
作为一个例子，考虑下面1001元素信号n x ，它包含了2个正弦信号和噪声平均功率通过用下述求和去近似积分求得你还可以用单边PSD去计算平均功率周期图性能下面从四个角度讨论周期图法估计的性能：泄漏，分辨率，偏差和方差。
频谱泄漏考虑有限长信号[]L x n ，把它表示成无限长序列[]x n 乘以一个有限长矩形窗[]R w n 的乘积的形式经常很有用：[][][]L R x n x n w n =⋅因为时域的乘积等效于频域的卷积，所以上式的傅立叶变换是()()()/2/21s s f L R sf X f X W f d f ρρρ-=-⎰前文中导出的表达式()()2ˆL xxs X f P f f L=说明卷积对周期图有影响。
正弦数据的卷积影响最容易理解。
假设[]x n 是M 个复正弦的和[]1k Mj n k k x n A e ω==∑其频谱是()()1Ms k k k X f f A f f δ==-∑对一个有限长序列，就变成了()()()()/211/21s s f M ML s k k R k R k k k sf X f f A f W f d A W f f f δρρρ==-=--=-∑∑⎰所以在有限长信号的频谱中，Dirac 函数被替换成了形式为()R k W f f -的项，该项对应于矩形窗的中心在k f 的频率响应。
一个矩形窗的频率响应形状是一个sinc 信号，如下所示该图显示了一个主瓣和若干旁瓣，最大旁瓣大约在主瓣下方处。
这些旁瓣说明了频谱泄漏效应。
无限长信号的功率严格的集中在离散频率点k f 处，而有限长信号在离散频率点k f 附近有连续的功率。
因为矩形窗越短，它的频率响应对Dirac 冲击的近似性越差，所以数据越短它的频谱泄漏越明显。
考虑下面的100个采样的序列注意到频谱泄露只视数据长度而定。
周期图确实只对有限数据样本进行计算，但是这和频谱泄露无关。
分辨率分辨率指的是区分频谱特征的能力，是分析谱估计性能的关键概念。
要区分两个在频率上离得很近的正弦，要求两个频率差大于任何一个信号泄漏频谱的主瓣宽度。
主瓣宽度定义为主瓣上峰值功率一半的点间的距离（3dB 带宽）。
该宽度近似等于/s f L两个频率为1f 2f 的正弦信号，可分辨条件是上例中频率间隔10Hz，数据长度要大于100抽才能使得周期图中两个频率可分辨。
下图是只有67个数据长度的情况上述对分辨率的讨论都是在高信噪比的情况进行的，因此没有考虑噪声。
当信噪比低的时候，谱特征的分辨更难，而且周期图上会出现一些噪声的伪像，如下所示估计偏差周期图是对PSD 的有偏估计。
期望值可以是()()()2/22/21ss f L xx R s s f X f E P W f d f L f L ρρρ-⎧⎫⎪⎪=-⎨⎬⎪⎪⎩⎭⎰ 该式和频谱泄漏中的()L X f 式相似，除了这里的表达式用的是平均功率而不是幅度。
这暗示了周期图产生的估计对应于一个有泄漏的PSD 而非真正的PSD 。
注意()2R W f ρ-本质上是一个三角Bartlett 窗（事实是两个矩形脉冲的卷积是三角脉冲。
）这导致了最大旁瓣峰值比主瓣峰值低27dB ，大致是非平方矩形窗的2倍。
周期图估计是渐进无偏的。
这从早期的一个观察结果可以明显看出，随着记录数据趋于无穷大，矩形窗对频谱对Dirac 函数的近似也就越来越好。
然而在某些情况下，周期图法估计很差劲即使数据够长，这是因为周期图法的方差，如下所述。
周期图法的方差()()()()222sin 2/var 1sin 2/L s xx s s X f Lf f P f f L L f f ππ⎧⎫⎡⎤⎛⎫⎪⎪⎢⎥≈+ ⎪⎨⎬ ⎪⎢⎥⎝⎭⎪⎪⎣⎦⎩⎭L 趋于无穷大，方差也不趋于0。
用统计学术语讲，该估计不是无偏估计。
然而周期图在信噪比大的时候仍然是有用的谱估计器，特别是数据够长。
Modified Periodogram 修正周期图法在fft 前先加窗，平滑数据的边缘。
可以降低旁瓣的高度。
旁瓣是使用矩形窗产生的陡峭的剪切引入的寄生频率，对于非矩形窗，结束点衰减的平滑，所以引入较小的寄生频率。
但是，非矩形窗增宽了主瓣，因此降低了频谱分辨率。
函数periodogram 允许指定对数据加的窗，例如默认的矩形窗和Hamming 窗事实上加Hamming窗后信号的主瓣大约是矩形窗主瓣的2倍。
对固定长度信号，Hamming 窗能达到的谱估计分辨率大约是矩形窗分辨率的一半。
这种冲突可以在某种程度上被变化窗所解决，例如Kaiser窗。
非矩形窗会影响信号的功率，因为一些采样被削弱了。
为了解决这个问题函数periodogram将窗归一化，有平均单位功率。
这样的窗不影响信号的平均功率。
修正周期图法估计的PSD是其中U是窗归一化常数假如U保证估计是渐进无偏的。
Welch法包括：将数据序列划分为不同的段（可以有重叠），对每段进行改进周期图法估计，再平均。
用对象，或pwelch函数。
默认情况下数据划分为4段，50%重叠，应用Hamming窗。
取平均的目的是减小方差，重叠会引入冗余但是加Hamming窗可以部分消除这些冗余，因为窗给边缘数据的权重比较小。