MATLAB计算信号的幅度谱与相位谱

MATLAB计算信号的幅度谱与相位谱

引用

CSDN

1.

https://blog.csdn.net/weixin_45333185/article/details/143101062

信号的频谱分析包括信号的幅度谱和相位谱，其中幅度谱是组成信号的各频谱成分的幅度，相位谱是组成信号的各频谱成分的初始相位。本文将详细介绍使用MATLAB计算信号幅度谱和相位谱的方法，特别指出在计算信号相位谱时如何避免相位信息误差。

一、MATLAB代码实现

幅度谱的计算方法已在之前的系列文章中详细介绍过，这里重点分析相位谱的计算。相位谱的计算公式如下：

计算信号幅度谱和相位谱的MATLAB代码如下：

%% 生成时域信号
f1 = 10;                       % 余弦波1的频率，单位Hz
f2 = 20;                       % 余弦波2的频率，单位Hz
p1 = 30;                       % 余弦波1的初始相位，单位度
p2 = 45;                       % 余弦波2的初始相位，单位度
A1 = 0.5;                      % 余弦波的振幅
A2 = 1;                        % 余弦波的振幅
fs = 320;                     % 采样率                   
Ts = 1/fs;                     % 采样时间间隔    
t  =0:1/fs:2-1/fs;             % 信号时长2秒
s1 = A1*cos(2*pi*f1*t + p1*pi/180);  % 信号1
s2 = A2*cos(2*pi*f2*t + p2*pi/180);  % 信号2
x = s1 + s2;                         % 和信号
figure()
subplot(2,2,1);
plot(t,x);
grid on;
title('time domain signal');
xlabel('t/s');
ylabel('x(t)');
ylim([-1.5 2])
%% 频域分析：幅度谱与相位谱
N = 512;                           % FFT size
X = 1/N*fftshift(fft(x,N));        % N-point complex DFT
df = fs/N;                         % 频率分辨率
n = -N/2:N/2-1;                    % FFT结果序号
f = n*df;                          % x轴索引转换为频率值
subplot(2,2,2);
plot(f,abs(X));                    % 幅度谱
grid on;
title('Magnitude spectrum');
xlabel('f (Hz)');
ylabel('|X(k)|');
xlim([-100 100])
% 直接计算相位谱
phase = atan2(imag(X),real(X))*180/pi;    % 原始相位信息是被计算误差污染的
subplot(2,2,3);
plot(f,phase);                            % 相位谱
grid on;
title('Phase spectrum (noisy)');
xlabel('f(Hz)');
ylabel('angle(X(k))');
xlim([-100 100])

二、仿真结果分析

信号的时域波形图：

信号的双边幅度谱：

信号的相位谱：

【相位谱分析】

从这个相位谱可以发现，相位信息误差导致不能得到正确的相位谱。

这个相位谱充满了噪声，几乎不能得出有效的相位信息。相位谱之所以噪声很大，是因为相位谱是通过计算FFT结果的虚部与实部之比的反正切得到的。这样，即使是微小的浮点舍入误差也会放大结果，并以相位信息的形式错误地表现出来。为了理解这一点，可以打印出FFT结果的前几个样本，可以发现它们并不完全等于0（虽然理论上应该等于0，但实际上所得到的值是非常小的数值，大约在10^-16的量级）。计算反正切将导致不正确的结果。解决方案是定义一个容差阈值，所有低于阈值的计算相位都将被忽略。下面给出正确的相位谱计算方法。

三、正确的相位谱计算方法与仿真结果

1.MATLAB代码

正确的相位谱计算的MATLAB代码如下：

% 排除信号误差对相位谱的影响
X2 = X;
threshold = max(abs(X))/10000;            % 设置误差门限
X2(abs(X)<threshold) = 0;                 % 小于门限的值清零
phase = atan2(imag(X2),real(X2))*180/pi;  % 提取相位信息
subplot(2,2,4);
plot(f,phase);                            % 相位谱
grid on;
title('Phase spectrum');
xlabel('f(Hz)');
ylabel('angle(X(k))');
xlim([-100 100])

2.正确的相位谱

所得到的相位谱如下：

分析：

正确的相位谱指出在频率f = 10Hz处的30°相位偏移，和在频率f = 20Hz处的45°相位偏移。相位谱是反对称的：在f = −10Hz处φ = −30°，在f = −20Hz处φ = −45°，这个仿真结果完全符合实值信号的理论分析。