使用fft查找每个谐波的相位 [英] Finding the phase of each harmonics using fft

问题描述

我使用Matlab.

我有一个正弦信号:

X(放大器:220/频率:50)

我要加上3个谐波:

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

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

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

我将所有信号加在一起(例如X包含3个谐波)，结果信号称为: Xt

I sum all the signals together (like X containing 3 harmonics), the resulting signal is called : 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 (幅度和频率)！

What I want to do is : find the 3 harmonics signal (their amplitude, frequency and phase) starting for the summed signal Xt and knowing the fundamental signal X (amplitude and frequency) !

到目前为止，我已经能够使用fft检索谐波的频率和幅度，现在的问题是找到谐波的相位(在我们的示例中为30°， 50°和90°).

So far, I was able using fft, to retrieve the frequencies and the amplitudes of the harmonics, the problem now is finding the phases of the harmonics (in our case : 30°, 50° and 90°).

推荐答案

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

The FFT returns you an array consisting of complex numbers. To define phases of the frequency components you need to use angle() function for complex numbers. Don't forget: the phase of your harmonics has to be given in radians.

这是代码:

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)');

这将导致混乱(但是您可以很好地看到振幅):

It will result in such a mess (but you can see your amplitudes very well):

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

You can see a lot of phase components on the second plot. But if you eliminate all the frequencies which correspond to zero amplitudes, you will see your phases.

我们在这里:

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()，所以:

Now you can see the phases, but all of them are shifted to 90 degrees. Why? Because the FFT works with cos() instead of sin(), so:

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

更新

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

What if the parameters of some signal components are not integer numbers?

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

Let's add a new component x4:

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

使用提供的代码，您将获得以下图表:

Using the provided code you will get the following plot:

这不是我们期望得到的，不是吗?问题在于频率样本的分辨率.该代码用谐波近似信号，该频率以1 Hz采样.显然，在77.77 Hz这样的频率下工作还不够.

This is not really what we expected to get, isn't it? The problem is in the resolution of the frequency samples. The code approximates the signal with harmonics, which frequencies are sampled with 1 Hz. It is obviously not enough to work with frequencies like 77.77 Hz.

频率分辨率等于信号时间的倒数.在我们之前的示例中，信号的长度为1秒，这就是为什么频率采样为1/1s=1Hz的原因.因此，为了提高分辨率，您需要扩展已处理信号的时间窗口.为此，只需更正可用的t:

The frequency resolution is equal to the inversed value of the signal's time. In our previous example the signal's length was 1 second, that's why the frequency sampling was 1/1s=1Hz. So in order to increase the resolution, you need to expand the time window of the processed signal. To do so just correct the definition of the vaiable t:

frq_res = 0.01; %desired frequency resolution

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

它将产生以下光谱:

更新2

没关系，必须分析哪个频率范围.信号分量可能来自非常高的范围，在下一个示例中显示.假设信号看起来像这样:

It does not matter, which frequency range has to be analyzed. The signal components can be from a very high range, what is shown in the next example. Suppose the signal looks like this:

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

这是结果图:

相位已移至-90度，这已在前面进行了解释.

The phases are shifted to -90 degrees, what was explained earlier.

这是代码:

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)');

这篇关于使用fft查找每个谐波的相位的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！