我用的是matlab。
我有一个正弦信号:
X (amp:220/频率:50)
我添加了 3 个谐波:
x1 => (h2) amp:30/Freq:100/phase:30°
x2 => (h4) amp:10/Freq:200/phase:50°
x3 => (h6) amp:05/Freq:300/phase:90°
我将所有信号加在一起(如包含 3 个谐波的 X)，产生的信号称为: 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;
我想要做的是:找到开始于求和信号 Xt 并知道基本信号 X (幅度和频率)的 3 个谐波信号(它们的幅度、频率和相位)!
到目前为止，我能够使用
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));
更新
如果某些信号分量的参数不是整数怎么办？
让我们添加一个新组件
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
它将产生以下光谱:
更新 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)');