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 2929i Xk 与 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 整个频谱图是以 Nyquist 频率为对称轴的 并且可以明显识别出信号中含有两种频率成分 15Hz 和 40Hz 由此可以知道 FFT 变换数 据的对称性 因此用 FFT 对信号做谱分析 只需考察 0 Nyquist 频率范围内的福频特性 若没有给出采样频率和采样间隔 则分析通常对归一化频率 0 1 进行 另外 振幅的大小 与所用采样点数有关 采用 128 点和 1024 点的相同频率的振幅是有不同的表现值 但在 同一幅图中 40Hz 与 15Hz 振动幅值之比均为 4 1 与真实振幅 0 5 2 是一致的 为了 与真实振幅对应 需要将变换后结果乘以 2 除以 N 例 2 x 0 5 sin 2 pi 15 t 2 sin 2 pi 40 t fs 100Hz 绘制 1 数据个数 N 32 FFT 所用的采样点数 NFFT 32 2 N 32 NFFT 128 3 N 136 NFFT 128 4 N 136 NFFT 512 clf fs 100 采样频率 Ndata 32 数据长度 N 32 FFT 的数据长度 n 0 Ndata 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 求取振幅 f 0 N 1 fs N 真实频率 subplot 2 2 1 plot f 1 N 2 mag 1 N 2 2 N 绘出 Nyquist 频率之前的振幅 xlabel 频率 Hz ylabel 振幅 title Ndata 32 Nfft 32 grid on Ndata 32 数据个数 N 128 FFT 采用的数据长度 n 0 Ndata 1 t n fs 时间序列 x 0 5 sin 2 pi 15 t 2 sin 2 pi 40 t y fft x N mag abs y f 0 N 1 fs N 真实频率 subplot 2 2 2 plot f 1 N 2 mag 1 N 2 2 N 绘出 Nyquist 频率之前的振幅 xlabel 频率 Hz ylabel 振幅 title Ndata 32 Nfft 128 grid on Ndata 136 数据个数 N 128 FFT 采用的数据个数 n 0 Ndata 1 t n fs 时间序列 x 0 5 sin 2 pi 15 t 2 sin 2 pi 40 t y fft x N mag abs y f 0 N 1 fs N 真实频率 subplot 2 2 3 plot f 1 N 2 mag 1 N 2 2 N 绘出 Nyquist 频率之前的振幅 xlabel 频率 Hz ylabel 振幅 title Ndata 136 Nfft 128 grid on Ndata 136 数据个数 N 512 FFT 所用的数据个数 n 0 Ndata 1 t n fs 时间序列 x 0 5 sin 2 pi 15 t 2 sin 2 pi 40 t y fft x N mag abs y f 0 N 1 fs N 真实频率 subplot 2 2 4 plot f 1 N 2 mag 1 N 2 2 N 绘出 Nyquist 频率之前的振幅 xlabel 频率 Hz ylabel 振幅 title Ndata 136 Nfft 512 grid on 结论 1 当数据个数和 FFT 采用的数据个数均为 32 时 频率分辨率较低 但没有由于添零而 导致的其他频率成分 2 由于在时间域内信号加零 致使振幅谱中出现很多其他成分 这是加零造成的 其振 幅由于加了多个零而明显减小 3 FFT 程序将数据截断 这时分辨率较高 4 也是在数据的末尾补零 但由于含有信号的数据个数足够多 FFT 振幅谱也基本不受 影响 对信号进行频谱分析时 数据样本应有足够的长度 一般 FFT 程序中所用数据点数与 原含有信号数据点数相同 这样的频谱图具有较高的质量 可减小因补零或截断而产生的 影响 例 3 x cos 2 pi 0 24 n cos 2 pi 0 26 n 1 数据点过少 几乎无法看出有关信号频谱的详细信息 2 中间的图是将 x n 补 90 个零 幅度频谱的数据相当密 称为高密度频谱图 但从图 中很难看出信号的频谱成分 3 信号的有效数据很长 可以清楚地看出信号的频率成分 一个是 0 24Hz 一个是 0 26Hz 称为高分辨率频谱 可见 采样数据过少 运用 FFT 变换不能分辨出其中的频率成分 添加零后可增加 频谱中的数据个数 谱的密度增高了 但仍不能分辨其中的频率成分 即谱的分辨率没有 提高 只有数据点数足够多时才能分辨其中的频率成分 第二 相谱 相位谱和频率普是回事儿 想着把频谱中的幅值部分换成相角就可以了 由于没有找到具体的理论 就举几个例子说明一下 比如要求 y sin 2 pi 60 t 的相位谱 程序如下 fs 200 N 1024 n 0 N 1 t n fs y sin 2 pi 60 t Y fft y N A abs Y f n fs N ph 2 angle Y 1 N 2 ph ph 180 pi plot f 1 N 2 ph 1 N 2 xlabel 频率 hz ylabel 相角 title 相位谱 grid on 期中的 ph 2 angle Y 1 N 2 ph ph 180 pi 是利用 angle 函数求出每个点的角度 并由 弧度转化成角度 angle 函数解释 Phase angle Syntax P angle Z Description P angle Z returns the phase angles in radians for each element of complex array Z The angles lie between For complex Z the magnitude R and phase angle theta are given by R abs Z theta angle Z and the statement Z R exp i theta converts back to the original complex Z Examples Z 1 1i 2 1i 3 1i 4 1i 1 2i 2 2i 3 2i 4 2i 1 3i 2 3i 3 3i 4 3i P angle Z P 0 7854 0 4636 0 3218 0 2450 1 1071 0 7854 0 5880 0 4636 1 2490 0 9828 0 7854 0 6435 1 3258 1 1071 0 9273 0 7854 Algorithms The angle function can be expressed as angle z imag log z atan2 imag z real z 第三 功率谱 matlab 实现经实现经典功率典功率谱谱估估计计 fft 做出来是频谱 psd 做出来是功率谱 功率谱丢失了频谱的相位信息 频谱不同的信号 其功率谱是可能相同的 功率谱是幅度取模后平方 结果是个实数 matlab 中自功率谱密度直接用 psd 函数就可以求 按照 matlab 的说法 psd 能实现 Welch 法估计 即相当于用改进的平均周期图法来求取随机信号的功率谱密度估计 psd 求出的 结果应该更光滑吧 1 直接法 直接法 直接法又称周期图法 它是把随机序列 x n 的 N 个观测数据视为一能量有限的序列 直接 计算 x n 的离散傅立叶变换 得 X k 然后再取其幅值的平方 并除以 N 作为序列 x n 真实功率谱的估计 Matlab 代码示例 clear Fs 1000 采样频率 n 0 1 Fs 1 产生含有噪声的序列 xn cos 2 pi 40 n 3 cos 2 pi 100 n randn size n window boxcar length xn 矩形窗 nfft 1024 Pxx f periodogram xn window nfft Fs 直接法 plot f 10 log10 Pxx 2 间间接法 接法 间接法先由序列 x n 估计出自相关函数 R n 然后对 R n 进行傅立叶变换 便得到 x n 的 功率谱估计 Matlab 代码示例 clear Fs 1000 采样频率 n 0 1 Fs 1 产生含有噪声的序列 xn cos 2 pi 40 n 3 cos 2 pi 100 n randn size n nfft 1024 cxn xcorr xn unbiased 计算序列的自相关函数 CXk fft cxn nfft Pxx abs CXk index 0 round nfft 2 1 k index Fs nfft plot Pxx 10 log10 Pxx index 1 plot k plot Pxx 3 改 改进进的直接法 的直接法 对于直接法的功率谱估计 当数据长度 N 太大时 谱曲线起伏加剧 若 N 太小 谱的分辨率又不好 因此需要改进 3 1 Bartlett 法法 Bartlett 平均周期图的方法是将 N 点的有限长序列 x n 分段求周期图再平均 Matlab 代码示例 clear Fs 1000 n 0 1 Fs 1 xn cos 2 pi 40 n 3 cos 2 pi 100 n randn size n nfft 1024 window boxcar length n 矩形窗 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 plot Pxxc 10 log10 Pxxc index 1 figure 1 plot k plot Pxx pause figure 2 plot k plot Pxx plot Pxx plot Pxxc plot Pxx plot Pxxc 3 2 Welch 法法 Welch 法对 Bartlett 法进行了两方面的修正 一是选择适当的窗函数 w n 并再周期图计 算前直接加进去 加窗的优点是无论什么样的窗函数均可使谱估计非负 二是在分段时 可使各段之间有重叠 这样会使方差减小 Matlab 代码示例 clear Fs 1000 n 0 1 Fs 1 xn cos 2 pi 40 n 3 cos 2 pi 100 n randn size n nfft 1024 window boxcar 100 矩形窗 window1 hamming 100 海明窗 window2 blackman 100 blackman 窗 noverlap 20 数据无重叠 range half 频率间隔为 0 Fs 2 只计算一半的频率 Pxx f pwelch xn window noverlap nfft Fs range Pxx1 f pwelch xn window1 noverlap nfft Fs range Pxx2 f pwelch xn window2 noverlap nfft Fs range plot Pxx 10 log10 Pxx plot Pxx1 10 log10 Pxx1 plot Pxx2 10 log10 Pxx2 figure 1 plot f plot Pxx pause figure 2 plot f plot Pxx1 pause figure 3 plot f plot Pxx2 第四 相关性分析 1 首先说说自相关和互相关的概念 这个是信号分析里的概念 他们分别表示的是两个时间序列之间和同一个时间序列在任 意两个不同时刻的取值之间的相关程度 即互相关函数是描述随机信号 x t y t 在任意两 个不同时刻 t1 t2 的取值之间的相关程度 自相关函数是描述随机信号 x t 在任意两个不 同时刻 t1 t2 的取值之间的相关程度 自相关函数是描述随机信号 X t 在任意两个不同时刻 t1 t2 的取值之间的相关程度 互相关函数给出了在频域内两个信号是否相关的一个 判断指标 把两测点之间信号的互谱与各自的自谱联系了起来 它能用来确定输出信号有 多大程度来自输入信号 对修正测量中接入噪声源而产生 的误差非常有效 事实上 在图象处理中 自相关和互相关函数的定义如下 设原函数是 f t 则自相关 函数定义为 R u f t f t 其中 表示卷积 设两个函数分别是 f t 和 g t 则互相关函 数定义为 R u f t g t 它反映的是两个函数在不同的相对位置上互相匹配的程度 那么 如何在 matlab 中实现这两个相关并用图像显示出来呢 dt 1 t 0 dt 100 x cos t a b xcorr x unbiased plot b dt a 上面代码是求自相关函数并作图 对于互相关函数 稍微修改一下就可以了 即把 a b xcorr x unbiased 改为 a b xcorr x y unbiased 便可 2 实现过程 在 Matalb 中 求解 xcorr 的过程事实上是利用 Fourier 变换中的卷积定理进行的 即 R u ifft fft f fft g 其中 表示乘法 注 此公式仅表示形式计算 并非实际计算所用 的公式 当然也可以直接采用卷积进行计算 但是结果会与 xcorr 的不同 事实上 两者 既然有定理保证 那么结果一定是相同的 只是没有用对公式而已 下面是检验两者结果 相同的代码 dt 1 t 0 dt 100 x 3 sin t y cos 3 t subplot 3 1 1 plot t x subplot 3 1 2 plot t y a b xcorr x y subplot 3 1 3 plot b dt a yy cos 3 fliplr t or use yy fliplr y z conv x yy pause subplot 3 1 3 plot b dt z r 即在 xcorr 中不使用 scaling 3 其他相关问题 1 相关程度与相关函数的取值有什么联系 相关系数只是一个比率 不是等单位量度 无什么单位名称 也不是相关的百分数 一 般取小数点后两位来表示 相关系数的正负号只表 示相关的方向 绝对值表示相关的程度 因为不是等单位的度量 因而不能说相关系数 0 7 是 0 35 两倍 只能说相关系数为 0 7 的二列变量相关程度 比相关系数为 0 35 的二列变量相关程度更为密切和更高 也不能说相关系数从 0 70 到 0 80 与相关系数从 0 30 到 0 40 增加的程度一样