Finding the phase of each harmonics using fft

4
votes

I use Matlab.

I have a sinusoidal signal :

X (amp:220/ Freq:50)

to which I add 3 harmonics :

x1 => (h2) amp:30 / Freq:100 / phase:30°

x2 => (h4) amp:10 / Freq:200 / phase:50°

x3 => (h6) amp:05 / Freq:300 / phase:90°

I sum all the signals together (like X containing 3 harmonics), the resulting signal is called : Xt

Here is the code :

%% 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;

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) !

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°).

2
matlabfftphase
Oh for the record, i'm new to matlab and i'm thinking : x1 = 30.*sin(2 .* pi .* 100 .* t + 30); is wrong because i can't use degree in matlab but radians ?! an i wrong ?! – KADEM Mohammed
Look at the 'arg' of the complex fft results. – greggo
i know, i used the abs to find the Amplitudes but i don't know how to use the the arg for find the phase ! – KADEM Mohammed
The value you get depends on how the signal is positioned in the fft input. Try adjusting one of the phases a little and see what happens. And the input to sin/cos, and output from arg, are in radians. – greggo
Phase relative to what point in time? Phase continuously changes unless you establish a fixed reference point for analysis, and that phase will be different unless the reference point or time is the same as the one used for synthesis or generation. – hotpaw2

2 Answers

6
votes

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.

Here is the code:

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):

enter image description here

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.

Here we are:

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

enter image description here

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

UPDATE

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

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:

fft of a signal with real number parameters

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.

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

It will result in the following spectra:

spectra for the signals with non-integer parameters

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

Here is the resulting plot:

enter image description here

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

Here is the code:

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

To start off we should note (as you correctly found out in comments) that Matlab uses radians for angles, so the harmonics should be:

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

The simple case

The process of estimating the amplitude, frequency and phase of the frequency components will usually start off with taking the Fast Fourier Transform (FFT), and selecting the strongest frequency components:

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

Now if the fundamental frequency and the frequency of the harmonics happens to be exact multiples of fs/N where fs is the sampling rate and N is the number of samples (in this case length(Xt)) then the tones will fall exactly on a bin, and the frequencies, amplitudes and phases of each component can be estimated fairly easily with:

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

The usual and more complicated reality

If on the other hand the frequency components are not exact multiples of fs/N, (or at least are not known to be exact multiples of fs/N, you are after all trying to estimate the frequency of those components) then things get more complicated. Note that this can have a particularly significant effect on the phase estimate.

We start off by recalling that a pure complex tone (exp(2*pi*j*n*f/fs)) of finite length N has a Discrete Fourier Transform (DFT) given by:

enter image description here

One estimation approach could be to start by estimating the frequency. The amplitude and phase can be factored out by looking at the ratio of the magnitudes of two successive bins of Xf around the peak, mainly at indices idx(i) and idx(i)+1. Under the assumption that those two bins suffer little interference, then the ratio can be expressed as:

ratio = abs(Xf(idx(i)+1)/Xf(idx)) 
      = abs(sin(pi*frac/N)/sin(pi*(frac-1)/N))

Where the frequency to be estimated is f = (idx(i)-1 + frac)*fs/N. The parameter frac can then be obtained with the Newton-Raphson method:

% 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

Which we use to estimate the frequency with:

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

Now that we have the tone frequency, we can obtain the amplitude and phase by fitting the DFT equation given above (where we also add a factor of 2 for the amplitude by since we are dealing with a real tone):

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

And putting it all together:

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