使用 fft 查找每个谐波的相位
Finding the phase of each harmonics using fft
我用的是Matlab
我有一个正弦信号：
X (amp:220/ Freq: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°
我将所有信号加在一起（比如 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;
我想要做的是：找到从求和信号 Xt 开始的 3 次谐波信号（它们的幅度、频率和相位）并知道基波信号 X（振幅和频率）!
到目前为止，我能够使用 fft 检索谐波的频率和振幅，现在的问题是找到谐波的相位（在我们的例子中：30°、50° 和 90°）。
FFT returns 你是一个由复数组成的数组。要定义频率分量的相位，您需要对复数使用 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 的原因。所以为了提高分辨率，需要扩大处理信号的时间window。为此，只需更正变量 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)');
首先我们应该注意（正如您在评论中正确发现的那样）Matlab 使用弧度表示角度，因此谐波应该是：
%% Harmonics
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);
简单案例
估计频率分量的幅度、频率和相位的过程通常从采用快速傅立叶变换 (FFT) 开始，然后选择最强的频率分量：
% Compute the frequency spectrum
N = length(Xt);
Xf = fft(Xt);
Nmax = N/2 + 1;
Xf = Xf(1:Nmax);
% Locate the peaks
largest_peak = max(20*log10(abs(Xf)));
peak_floor = largest_peak - 100; % to reject peaks from spectral leakage and noise
[pks,idx] = findpeaks((max(peak_floor, 20*log10(abs(Xf))) - peak_floor)')
现在如果基频和谐波频率正好是fs/N的倍数，其中fs是采样率，N是样本数（在在这种情况下 length(Xt)) 那么音调将恰好落在一个 bin 上，并且每个分量的频率、幅度和相位都可以很容易地估计为：
Amp = 2*abs(Xf(idx))/N;
freq = (idx-1)*fs/N;
phase = angle(Xf(idx));
phase = phase - phase(1); % set phase reference to that of the fundamental
平凡而复杂的现实
另一方面，如果频率分量不是 fs/N 的精确倍数，（或者至少不知道是 fs/N 的精确倍数，你毕竟是在尝试估计这些组件的频率）然后事情变得更加复杂。请注意，这会对相位估计产生特别显着的影响。
我们首先回忆一下有限长度 N 的纯复音 (exp(2*pi*j*n*f/fs)) 的离散傅里叶变换 (DFT) 由下式给出：
一种估计方法可以从估计频率开始。幅度和相位可以通过查看峰值周围 Xf 的两个连续 bin 的幅度之比来分解，主要在索引 idx(i) 和 idx(i)+1 处。假设这两个 bin 受到的干扰很小，则该比率可以表示为：
ratio = abs(Xf(idx(i)+1)/Xf(idx))
= abs(sin(pi*frac/N)/sin(pi*(frac-1)/N))
其中要估计的频率是 f = (idx(i)-1 + frac)*fs/N。然后可以用Newton-Raphson方法得到参数frac：
% Solve for "f" for which ratio = sin(pi*frac/N)/sin(pi*(frac-1)/N)
function f = fractional_frequency(ratio, N)
niter = 20;
K = (pi/N) * sin(pi/N);
f = 0;
for i=1:niter
a = sin(pi*f/N);
b = sin(pi*(f-1)/N);
y = ratio - a/b;
yp = K / (b^2);
f = max(-0.5, min(f - y/yp, 0.5));
end
end
我们用它来估计频率：
freq = zeros(1,length(idx));
for i=1:length(idx)
ratio = abs(Xf(idx(i)+1))/abs(Xf(idx(i)));
if (abs(Xf(idx(i)+1)) > abs(Xf(idx(i)-1)))
ratio = -ratio;
end
frac = fractional_frequency(ratio, N)
freq(i) = (idx(i)-1+frac)*fs/N;
end
现在我们有了音调频率，我们可以通过拟合上面给出的 DFT 方程来获得振幅和相位（由于我们处理的是真实音调，我们还为振幅添加了一个因子 2）：
Amp(i) = 2 * abs(Xf(idx(i))) * abs(sin(pi*frac/N)/sin(pi*frac));
phase(i) = angle( Xf(idx(i)) .* (1-exp(2*pi*frac*j/N)) ./ (1-exp(2*pi*frac*j)) );
把它们放在一起：
Amp = zeros(1,length(idx));
freq = zeros(1,length(idx));
phase = zeros(1,length(idx));
for i=1:length(idx)
ratio = abs(Xf(idx(i)+1))/abs(Xf(idx(i)));
if (abs(Xf(idx(i)+1)) > abs(Xf(idx(i)-1)))
ratio = -ratio;
end
frac = fractional_frequency(ratio, N)
freq(i) = (idx(i)-1+frac)*fs/N;
Amp(i) = 2 * abs(Xf(idx(i))) * abs(sin(pi*frac/N)/sin(pi*frac));
phase(i) = angle( Xf(idx(i)) .* (1-exp(2*pi*frac*j/N)) ./ (1-exp(2*pi*frac*j)) );
end
phase = phase - phase(1); % set phase reference to that of the fundamental
我用的是Matlab
我有一个正弦信号：
X (amp:220/ Freq: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°
我将所有信号加在一起（比如 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;
我想要做的是：找到从求和信号 Xt 开始的 3 次谐波信号（它们的幅度、频率和相位）并知道基波信号 X（振幅和频率）!
到目前为止，我能够使用 fft 检索谐波的频率和振幅，现在的问题是找到谐波的相位（在我们的例子中：30°、50° 和 90°）。
FFT returns 你是一个由复数组成的数组。要定义频率分量的相位，您需要对复数使用 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 的原因。所以为了提高分辨率，需要扩大处理信号的时间window。为此，只需更正变量 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)');
首先我们应该注意（正如您在评论中正确发现的那样）Matlab 使用弧度表示角度，因此谐波应该是：
%% Harmonics
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);
简单案例
估计频率分量的幅度、频率和相位的过程通常从采用快速傅立叶变换 (FFT) 开始，然后选择最强的频率分量：
% Compute the frequency spectrum
N = length(Xt);
Xf = fft(Xt);
Nmax = N/2 + 1;
Xf = Xf(1:Nmax);
% Locate the peaks
largest_peak = max(20*log10(abs(Xf)));
peak_floor = largest_peak - 100; % to reject peaks from spectral leakage and noise
[pks,idx] = findpeaks((max(peak_floor, 20*log10(abs(Xf))) - peak_floor)')
现在如果基频和谐波频率正好是fs/N的倍数，其中fs是采样率，N是样本数（在在这种情况下 length(Xt)) 那么音调将恰好落在一个 bin 上，并且每个分量的频率、幅度和相位都可以很容易地估计为：
Amp = 2*abs(Xf(idx))/N;
freq = (idx-1)*fs/N;
phase = angle(Xf(idx));
phase = phase - phase(1); % set phase reference to that of the fundamental
平凡而复杂的现实
另一方面，如果频率分量不是 fs/N 的精确倍数，（或者至少不知道是 fs/N 的精确倍数，你毕竟是在尝试估计这些组件的频率）然后事情变得更加复杂。请注意，这会对相位估计产生特别显着的影响。
我们首先回忆一下有限长度 N 的纯复音 (exp(2*pi*j*n*f/fs)) 的离散傅里叶变换 (DFT) 由下式给出：
一种估计方法可以从估计频率开始。幅度和相位可以通过查看峰值周围 Xf 的两个连续 bin 的幅度之比来分解，主要在索引 idx(i) 和 idx(i)+1 处。假设这两个 bin 受到的干扰很小，则该比率可以表示为：
ratio = abs(Xf(idx(i)+1)/Xf(idx))
= abs(sin(pi*frac/N)/sin(pi*(frac-1)/N))
其中要估计的频率是 f = (idx(i)-1 + frac)*fs/N。然后可以用Newton-Raphson方法得到参数frac：
% Solve for "f" for which ratio = sin(pi*frac/N)/sin(pi*(frac-1)/N)
function f = fractional_frequency(ratio, N)
niter = 20;
K = (pi/N) * sin(pi/N);
f = 0;
for i=1:niter
a = sin(pi*f/N);
b = sin(pi*(f-1)/N);
y = ratio - a/b;
yp = K / (b^2);
f = max(-0.5, min(f - y/yp, 0.5));
end
end
我们用它来估计频率：
freq = zeros(1,length(idx));
for i=1:length(idx)
ratio = abs(Xf(idx(i)+1))/abs(Xf(idx(i)));
if (abs(Xf(idx(i)+1)) > abs(Xf(idx(i)-1)))
ratio = -ratio;
end
frac = fractional_frequency(ratio, N)
freq(i) = (idx(i)-1+frac)*fs/N;
end
现在我们有了音调频率，我们可以通过拟合上面给出的 DFT 方程来获得振幅和相位（由于我们处理的是真实音调，我们还为振幅添加了一个因子 2）：
Amp(i) = 2 * abs(Xf(idx(i))) * abs(sin(pi*frac/N)/sin(pi*frac));
phase(i) = angle( Xf(idx(i)) .* (1-exp(2*pi*frac*j/N)) ./ (1-exp(2*pi*frac*j)) );
把它们放在一起：
Amp = zeros(1,length(idx));
freq = zeros(1,length(idx));
phase = zeros(1,length(idx));
for i=1:length(idx)
ratio = abs(Xf(idx(i)+1))/abs(Xf(idx(i)));
if (abs(Xf(idx(i)+1)) > abs(Xf(idx(i)-1)))
ratio = -ratio;
end
frac = fractional_frequency(ratio, N)
freq(i) = (idx(i)-1+frac)*fs/N;
Amp(i) = 2 * abs(Xf(idx(i))) * abs(sin(pi*frac/N)/sin(pi*frac));
phase(i) = angle( Xf(idx(i)) .* (1-exp(2*pi*frac*j/N)) ./ (1-exp(2*pi*frac*j)) );
end
phase = phase - phase(1); % set phase reference to that of the fundamental