Phase extraction from FFT

Hello everyone,

I have written a code that takes as input voltage and current waveforms and is intended to have as output their frequency components and the power consumed. The thing is, by having a really high sampling frequency, tha data files must be very big. So I used zero-padding to fill in for the missing data. My comment is that I am getting wrong phase results and some of the magnitudes differ slightly from the real ones. What is really troubling though is the wrong phase results.

Below is my code:.

1. This is a simplified version of the code which is intended to help me find my mistakes. The actual code gets input from data sheets.

Could you please inform me if something is wrong with my code or if I am missing something?

clear all;
clc;
F=13.56e6;
t=-0.5e-5:1/(2.5e9):0.5e-5;
V=-500+230*sin(2*pi*F*t-pi/4)+1*sin(2*pi*2*F*t-pi/6)+0.5*sin(2*pi*3*F*t-pi/3);
I=0.1*sin(2*pi*F*t+pi/8)+0.02*sin(2*pi*2*F*t+pi/10)+0.04*sin(2*pi*3*F*t);

signal1=V;
signal2=I;
N=length(t);  
dt=-t(1,1)+t(1,2);

NFFT=5e5;
FS=1/dt; 
Fn=FS/2;                         %nyquist frequency

Y1=fft(signal1,NFFT);
Y2=fft(signal2,NFFT);

F=linspace(0,1,(fix(NFFT/2)+1))*Fn;    %freq vector

Iv=1:length(F);

HAmpl1=2*abs(Y1(Iv))/N;
HAmpl1(1,1)=HAmpl1(1,1)/2;
HAmpl2=2*abs(Y2(Iv))/N;
phase1= rad2deg(angle(Y1(Iv)));
phase2= rad2deg(angle(Y2(Iv)));

%% Figuring where the harmonics are in the matrix

C=zeros(1,length(F));
YY=zeros(1,length(Y1));
YYY=zeros(1,length(Y1));
Num=3;                  %number of harmonics calculated
M=max(HAmpl1(1,max(find(F<1.2e7 & F>1e7)):end));

noharm=ones(1,Num+1);
noharm(2)=find(HAmpl1==M);
temp=noharm(2);
for i=2:Num
    noharm(i+1)=i*temp-(i-1);
end

%Values of interest
HAmplinterestV(1,:)=HAmpl1(noharm);
PhaseinterestV(1,:)=phase1(noharm);
HAmplinterestC(1,:)=HAmpl2(noharm);
PhaseinterestC(1,:)=phase2(noharm);

figure(1)
Ampl=subplot(2,1,1)
stem(F(noharm),HAmpl1(noharm))
Phasee=subplot(2,1,2)
stem(F(noharm),phase1(noharm))

[Rafael Hernandez-Walls]

What values ​​does the variable tt take?

[Bill G.]

Rafael Hernandez-Walls My apologies, it was left over from another part of the code, but it shouldn't matter. I was recreating the signal from the harmonic content by selecting certain harmonic frequencies and then impelementing ifft, my way of verifying the results.

[Bill G.]

Star Strider Well, using unwrap will render my intuitive error recognition useless, at least if I understood it's purpose correctly. Nevertheless, shouldn't the phases of let's say the voltage, be close to or even equal to pi/4(45 deg) for the fundamental, pi/6(30 deg) for the 2nd harmonic and pi/3(60 deg) for the 3rd?

[Star Strider]

The waveform is (at the very least) not ‘pure’. Waveforms with multiple frequencies have multiple phase relationships. Most power calculations (that I am aware of) involve single frequencies. Dealing with waveform such as these are analytically difficult. Even though they involve relatively elementary trigonometry, interpreting the results can be a challenge.

% clear all;
% clc;
F=13.56e6;
t=-0.5e-5:1/(2.5e9):0.5e-5;
V=-500+230*sin(2*pi*F*t-pi/4)+1*sin(2*pi*2*F*t-pi/6)+0.5*sin(2*pi*3*F*t-pi/3);
I=0.1*sin(2*pi*F*t+pi/8)+0.02*sin(2*pi*2*F*t+pi/10)+0.04*sin(2*pi*3*F*t);

signal1=V;
signal2=I;
N=length(t);  
dt=-t(1,1)+t(1,2);

NFFT=5e5;
FS=1/dt; 
Fn=FS/2;                         %nyquist frequency

Y1=fft(signal1,NFFT);
Y2=fft(signal2,NFFT);

F=linspace(0,1,(fix(NFFT/2)+1))*Fn;    %freq vector

Iv=1:length(F);

HAmpl1=2*abs(Y1(Iv))/N;
HAmpl1(1,1)=HAmpl1(1,1)/2;
HAmpl2=2*abs(Y2(Iv))/N;
phase1= rad2deg(angle(Y1(Iv)));
phase2= rad2deg(angle(Y2(Iv)));

%% Figuring where the harmonics are in the matrix

C=zeros(1,length(F));
YY=zeros(1,length(Y1));
YYY=zeros(1,length(Y1));
Num=3;                  %number of harmonics calculated
M=max(HAmpl1(1,max(find(F<1.2e7 & F>1e7)):end));

noharm=ones(1,Num+1);
noharm(2)=find(HAmpl1==M);
temp=noharm(2);
for i=2:Num
    noharm(i+1)=i*temp-(i-1);
end

%Values of interest
HAmplinterestV(1,:)=HAmpl1(noharm);
PhaseinterestV(1,:)=phase1(noharm);
HAmplinterestC(1,:)=HAmpl2(noharm);
PhaseinterestC(1,:)=phase2(noharm);

figure(1)
Ampl=subplot(2,1,1)

Ampl = 
  Axes with properties:

             XLim: [0 1]
             YLim: [0 1]
           XScale: 'linear'
           YScale: 'linear'
    GridLineStyle: '-'
         Position: [0.1300 0.5838 0.7750 0.3412]
            Units: 'normalized'

  Show all properties

stem(F(noharm),HAmpl1(noharm))
Phasee=subplot(2,1,2)

Phasee = 
  Axes with properties:

             XLim: [0 1]
             YLim: [0 1]
           XScale: 'linear'
           YScale: 'linear'
    GridLineStyle: '-'
         Position: [0.1300 0.1100 0.7750 0.3412]
            Units: 'normalized'

  Show all properties

stem(F(noharm),phase1(noharm))

figure
yyaxis left
plot(t,V)
ylabel('V(t)')
yyaxis right
plot(t,I)
ylabel('I(t)')
grid
xlim([-1 1]*2.5E-7)

figure
plot(t, V.*I)
grid
title('Power')
xlim([-1 1]*2.5E-7)

.

[Bill G.]

Thank you very much for the reply.

So, is there a way to get the amplitude and phase of each harmonic as i put it in the formula above? For example,

I put 230*sin(2*pi*F*t-pi/4). Can I get as phase ~pi/4 and amplitude ~230 somehow?

Oh and btw if I make the time vector to be of length 500k points, which is the dimension of my FFT, phases are projected to their exact value. Is there any way to make them correct with data size of 25001 points (what I am using above).