Thread Subject:

reconstruct waveformm from harmonic of FFT

I have a noisy signal assumed to be a sum of odd harmonics. I want to use the amplitudes and phase angles from the FFT to reconstruct a smooth version.

Here is my code:

clear all
fs=50; % sampling frequency

t=0:1/fs:12345/fs; % the length of signal is 12345.
t=t';

x1=sin(2*pi*t); % the input signal is a sinusoidal wave;
x2=12*sin(2*pi*t+1)+3*sin(3*2*pi*t+3)+0.8*sin(5*2*pi*t+2)+0.02*sin(7*2*pi*t+1)+0.05*randn(size(t)); % the output is a sum of odd harmonics

x=[x1 x2];


fi=sinfapm(x(:,1),fs); % the frequency is obtained by sinewave fitting.

Nr_cycles=floor(fi*length(x)/fs); % number of comlete cycles

N=round(Nr_cycles*fs/fi); % cutoff length for coherent sampling
x=x(1:N,:);

t=t(1:N); % cut the time axis too

Xraw=fft(x(:,2));
phase=angle(Xraw);

freq_bin=0:N-1; % frequency bin
freq_res=fs/N; % frequency resolution
freq_i=freq_bin*freq_res; % frequency axis

cutoff=ceil(N/2); % cut the first half of FFT
Xraw=Xraw(1:cutoff);
phase=phase(1:cutoff);


max_nr_harm=7; % in practice the max number of harmonics is determined by the Nyquist frequency

I_n=zeros((max_nr_harm+1)/2,1); % amplitudes of the odd harmonics
delta_n=zeros((max_nr_harm+1)/2,1); % phase angles of the odd harmonics

% extract the harmonic information
for n=1:(max_nr_harm+1)/2
    I_n(n)=abs(Xraw(round((2*n-1)*fi/freq_res)+1));
    delta_n(n)=phase(round((2*n-1)*fi/freq_res)+1);
end

reconst=zeros(length(x),1);

% reconstruct the signal by summing up the sine components
for n=1:2:max_nr_harm
    reconst=reconst+I_n((n+1)/2)*sin(n*2*pi*fi*t+delta_n((n+1)/2));
end

plot(t,x(:,2),'k',t,reconst,'r')


The constructed result is totally wrong compared with the original one. I noticed that the amplitudes of the odd harmonics are correctly extracted, but the phase angles are wrong. If I substitute the correct phase angles the reconstructed signal is perfect. How to get the phase angles correct?

Thank you in advance!

"Andrew Sun" <andrewx100@gmail.com> wrote in message <jpthbr$7qa$1@newscl01ah.mathworks.com>...
> I have a noisy signal assumed to be a sum of odd harmonics. I want to use the amplitudes and phase angles from the FFT to reconstruct a smooth version.
>
> Here is my code:
>
> clear all
> fs=50; % sampling frequency
>
> t=0:1/fs:12345/fs; % the length of signal is 12345.
> t=t';
>
> x1=sin(2*pi*t); % the input signal is a sinusoidal wave;
> x2=12*sin(2*pi*t+1)+3*sin(3*2*pi*t+3)+0.8*sin(5*2*pi*t+2)+0.02*sin(7*2*pi*t+1)+0.05*randn(size(t)); % the output is a sum of odd harmonics
>
> x=[x1 x2];
>
>
> fi=sinfapm(x(:,1),fs); % the frequency is obtained by sinewave fitting.
>
> Nr_cycles=floor(fi*length(x)/fs); % number of comlete cycles
>
> N=round(Nr_cycles*fs/fi); % cutoff length for coherent sampling
> x=x(1:N,:);
>
> t=t(1:N); % cut the time axis too
>
> Xraw=fft(x(:,2));
> phase=angle(Xraw);
>
> freq_bin=0:N-1; % frequency bin
> freq_res=fs/N; % frequency resolution
> freq_i=freq_bin*freq_res; % frequency axis
>
> cutoff=ceil(N/2); % cut the first half of FFT
> Xraw=Xraw(1:cutoff);
> phase=phase(1:cutoff);
>
>
> max_nr_harm=7; % in practice the max number of harmonics is determined by the Nyquist frequency
>
> I_n=zeros((max_nr_harm+1)/2,1); % amplitudes of the odd harmonics
> delta_n=zeros((max_nr_harm+1)/2,1); % phase angles of the odd harmonics
>
> % extract the harmonic information
> for n=1:(max_nr_harm+1)/2
> I_n(n)=abs(Xraw(round((2*n-1)*fi/freq_res)+1));
> delta_n(n)=phase(round((2*n-1)*fi/freq_res)+1);
> end
>
> reconst=zeros(length(x),1);
>
> % reconstruct the signal by summing up the sine components
> for n=1:2:max_nr_harm
> reconst=reconst+I_n((n+1)/2)*sin(n*2*pi*fi*t+delta_n((n+1)/2));
> end
>
> plot(t,x(:,2),'k',t,reconst,'r')
>
>
> The constructed result is totally wrong compared with the original one. I noticed that the amplitudes of the odd harmonics are correctly extracted, but the phase angles are wrong. If I substitute the correct phase angles the reconstructed signal is perfect. How to get the phase angles correct?
>
> Thank you in advance!

Sorry! I forgot to scale the fft in the above code, which should be scale as:

Xraw=Xraw*2/N;

But this does not affect the problem.

"Andrew Sun" <andrewx100@gmail.com> wrote in message <jpthbr$7qa$1@newscl01ah.mathworks.com>...
>
>
How to get the phase angles correct?
==========

Add pi/2 to all of them.

"Matt J" wrote in message <jptp3b$743$1@newscl01ah.mathworks.com>...
> "Andrew Sun" <andrewx100@gmail.com> wrote in message <jpthbr$7qa$1@newscl01ah.mathworks.com>...
> >
> >
> How to get the phase angles correct?
> ==========
>
> Add pi/2 to all of them.

Problem solved. Thank you very much! Though I still wonder why there is a pi/2 lag...

"Andrew Sun" <andrewx100@gmail.com> wrote in message <jpuh43$7vg$1@newscl01ah.mathworks.com>...
>
> > Add pi/2 to all of them.
>
> Problem solved. Thank you very much! Though I still wonder why there is a pi/2 lag...
=============

Remember that an FFT doesn't decompose your signal into sine waves. It decomposes it into complex exponentials. The two are related according to the identity

sin(x)=exp(j*x)/(2j) - exp(-j*x)/(2j)

From this you can see that the complex exponentials with positive frequencies have phases pi/2 less than the sine wave (due to the 1/j coefficients).

"Matt J" wrote in message <jpvpc7$bu9$1@newscl01ah.mathworks.com>...
> "Andrew Sun" <andrewx100@gmail.com> wrote in message <jpuh43$7vg$1@newscl01ah.mathworks.com>...
> >
> > > Add pi/2 to all of them.
> >
> > Problem solved. Thank you very much! Though I still wonder why there is a pi/2 lag...
> =============
>
> Remember that an FFT doesn't decompose your signal into sine waves. It decomposes it into complex exponentials. The two are related according to the identity
>
> sin(x)=exp(j*x)/(2j) - exp(-j*x)/(2j)
>
> From this you can see that the complex exponentials with positive frequencies have phases pi/2 less than the sine wave (due to the 1/j coefficients).

I see. Thank you! I should have ponder more on the principle of DFT.