MATLAB: Extracting phase informations from a FFT

fftharmonicsphase

Hi,
I am currently struggling with a phase information extraction issue from a FFT. Here is basically what I want to do…
I have a signal coming from measurements that I want to replicate using the first 8 harmonics or so. I'm using the FFT to extract the amplitude of those harmonics which is pretty straight forward.
The problem occurs when I try to extract the phase of those harmonics from the same FFT data. I understood by reading some questions/answers here that it is not possible to extract directly the phase from a FFT except if some precises conditions are met(for instance, to have exactly one period of the signal without any overlapping, which is not the case here).
Is there a general way to get the information on the phase from the information I have in hand?
I already tried to inject a sinusoidal at a known frequency and phase in order to measure the phase difference between this signal and each of the harmonics instead of trying to extract the absolute phase, but it didn't work out well…
Any tips or any possible solution would be greatly appreciated!
Thank you, Alex

Best Answer

If you know the frequencies, you can formulate this problem as linear regression problem and obtain the least squares estimates.
1.) Form a design matrix where the first column is a column of ones (for the DC term-- or mean value). Then you need two columns for every harmonic, one cosine and one sine.
For example:
    t = linspace(0,1,1000)';
    x = cos(2*pi*100*t-pi/4)+1/2*cos(2*pi*200*t-pi/8)+ 0.1*randn(size(t));    
    X = ones(1000,5);
    X(:,2) = cos(2*pi*100*t)';
    X(:,3) = sin(2*pi*100*t)';
    X(:,4) = cos(2*pi*200*t)';
    X(:,5) = sin(2*pi*200*t)';
    beta = X\x;
    ampest100 = sqrt(beta(2)^2+beta(3)^2);
    phase100 = atan2(-beta(3),beta(2));
    ampest200 = sqrt(beta(4)^2+beta(5)^2);
    phase200 = atan2(-beta(5),beta(4));
Note the amplitude of the 100 Hz component is 1 and the phase is -pi/4. Note the amplitude of the 200 Hz component is 1/2 and the phase is -pi/8
Look how close ampest100 is to the amplitude of the 100 Hz component. Low how close phase100 is to the phase of the 100 Hz component.
Do the same for the 200 Hz component.
Related Question