使用 fft 查找每个谐波的相位

我用的是Matlab。

我有一个正弦信号：

X（放大器：220/频率：50）

我添加了 3 个谐波：

x1 => (h2) 放大器:30 / 频率:100 / 相位:30°

x2 => (h4) 放大器:10 / 频率:200 / 相位:50°

x3 => (h6) 放大器:05 / 频率:300 / 相位:90°

我将所有信号加在一起（例如 X 包含 3 个谐波），所得信号称为：Xt

这是代码：

%% Original signal
X = 220.*sin(2 .* pi .* 50 .* t);

%% Harmonics
x1 = 30.*sin(2 .* pi .* 100 .* t + 30);
x2 = 10.*sin(2 .* pi .* 200 .* t + 50);
x3 = 05.*sin(2 .* pi .* 300 .* t + 90);

%% adding the harmonics
Xt = X + x1 + x2 + x3;

我想要做的是：找到从求和信号开始的 3 个谐波信号（它们的幅度、频率和相位）Xt并了解基本信号X（幅度和频率）！

到目前为止，我能够使用fft，检索谐波的频率和幅度，现在的问题是找到谐波的相位（在我们的例子中：30°、50° 和 90°）。

FFT 返回一个由复数组成的数组。要定义频率分量的相位，您需要对复数使用 angle() 函数。不要忘记：谐波的相位必须以弧度表示。

这是代码：

Fs = 1000; % Sampling frequency                     

t=0 : 1/Fs : 1-1/Fs; %time

X = 220*sin(2 * pi * 50 * t);

x1 = 30*sin(2*pi*100*t + 30*(pi/180));
x2 = 10*sin(2*pi*200*t + 50*(pi/180));
x3 = 05*sin(2*pi*300*t + 90*(pi/180));

%% adding the harmonics
Xt = X + x1 + x2 + x3;

%Transformation
Y=fft(Xt); %FFT

df=Fs/length(Y); %frequency resolution

f=(0:1:length(Y)/2)*df; %frequency axis


subplot(2,1,1);
M=abs(Y)/length(Xt)*2; %amplitude spectrum
stem(f, M(1:length(f)), 'LineWidth', 0.5);
xlim([0 350]);
grid on;  

xlabel('Frequency (Hz)')
ylabel('Magnitude');

subplot(2,1,2);
P=angle(Y)*180/pi; %phase spectrum (in deg.)
stem(f, P(1:length(f)), 'LineWidth', 0.5);
xlim([0 350]);
grid on;

xlabel('Frequency (Hz)');
ylabel('Phase (degree)');

它会导致如此混乱（但你可以很好地看到你的振幅）：

您可以在第二张图上看到很多相位分量。但是，如果您消除与零幅度相对应的所有频率，您将看到您的相位。

我们到了：

Y=fft(Xt); %FFT

df=Fs/length(Y); %frequency resolution

f=(0:1:length(Y)/2)*df; %frequency axis

subplot(2,1,1);
M=abs(Y)/length(Xt)*2; %amplitude spectrum

M_rounded = int16(M(1:size(f, 2))); %Limit the frequency range
ind = find(M_rounded ~= 0);

stem(f(ind), M(ind), 'LineWidth', 0.5);
xlim([0 350]);
grid on;  

xlabel('Frequency (Hz)')
ylabel('Magnitude');

subplot(2,1,2);
P=angle(Y)*180/pi; %phase spectrum (in deg.)
stem(f(ind), P(ind), 'LineWidth', 0.5);
xlim([0 350]);
ylim([-100 100]);
grid on;

xlabel('Frequency (Hz)');
ylabel('Phase (degree)');

现在您可以看到相位，但所有相位都偏移到 90 度。为什么？因为 FFT 使用 cos() 而不是 sin()，所以：

X = 220*sin(2*pi*50*t + 0*(pi/180)) = 220*cos(2*pi*50*t - 90*(pi/180));

UPDATE

如果某些信号分量的参数不是整数怎么办？

让我们添加一个新组件x4:

x4 = 62.75*cos(2*pi*77.77*t + 57.62*(pi/180));

使用提供的代码，您将得到以下图：

这并不是我们真正期望得到的，不是吗？问题在于频率样本的分辨率。该代码用谐波来近似信号，其频率以 1 Hz 采样。显然，仅使用 77.77 Hz 这样的频率是不够的。

频率分辨率等于信号时间的倒数。在我们前面的例子中，信号的长度是 1 秒，这就是频率采样的原因1/1s=1Hz。因此，为了提高分辨率，需要扩大处理信号的时间窗口。为此，只需更正可用的定义t:

frq_res = 0.01; %desired frequency resolution

t=0 : 1/Fs : 1/frq_res-1/Fs; %time

它将产生以下光谱：

UPDATE 2

必须分析哪个频率范围并不重要。信号分量可能来自非常高的范围，如下一个示例所示。假设信号如下所示：

f=20e4; % 200 KHz
Xt = sin(2*pi*(f-55)*t + pi/7) + sin(2*pi*(f-200)*t-pi/7);

这是结果图：

相位偏移至 -90 度，如前所述。

这是代码：

Fs = 300e4; % Sampling frequency                     

frq_res = 0.1; %desired frequency resolution

t=0 : 1/Fs : 1/frq_res-1/Fs; %time

f=20e4;
Xt = sin(2*pi*(f-55)*t + pi/7) + sin(2*pi*(f-200)*t-pi/7);

Y=fft(Xt); %FFT

df=Fs/length(Y); %frequency resolution

f=(0:1:length(Y)/2)*df; %frequency axis

subplot(2,1,1);
M=abs(Y)/length(Xt)*2; %amplitude spectrum

M_rounded = int16(M(1:size(f, 2))); %Limit the frequency range
ind = find(M_rounded ~= 0);

stem(f(ind), M(ind), 'LineWidth', 0.5);
xlim([20e4-300 20e4]);
grid on;  

xlabel('Frequency (Hz)')
ylabel('Magnitude');

subplot(2,1,2);
P=angle(Y)*180/pi; %phase spectrum (in deg.)
stem(f(ind), P(ind), 'LineWidth', 0.5);
xlim([20e4-300 20e4]);
ylim([-180 180]);
grid on;

xlabel('Frequency (Hz)');
ylabel('Phase (degree)');