Нахождение фазы каждой гармоники с помощью БПФ

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 Гц. Очевидно, что недостаточно работать с такими частотами, как 77,77 Гц. Разрешение по частоте равно обратному значению времени сигнала. В нашем предыдущем примере длина сигнала была равна 1 во-вторых, именно поэтому частотная выборка была 1/1s=1Hz. Поэтому для того, чтобы увеличить разрешение, необходимо расширить временное окно обрабатываемого сигнала. Для этого просто исправьте определение vaiable 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)');