MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Thread Subject: reconstruct waveformm from harmonic of FFT

 Subject: reconstruct waveformm from harmonic of FFT From: Andrew Sun Date: 27 May, 2012 15:31:07 Message: 1 of 6 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!
 Subject: reconstruct waveformm from harmonic of FFT From: Andrew Sun Date: 27 May, 2012 15:44:07 Message: 2 of 6 "Andrew Sun" wrote in message ... > 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.
 Subject: reconstruct waveformm from harmonic of FFT From: Matt J Date: 27 May, 2012 17:43:07 Message: 3 of 6 "Andrew Sun" wrote in message ... > > How to get the phase angles correct? ========== Add pi/2 to all of them.
 Subject: reconstruct waveformm from harmonic of FFT From: Andrew Sun Date: 28 May, 2012 00:33:07 Message: 4 of 6 "Matt J" wrote in message ... > "Andrew Sun" wrote in message ... > > > > > 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...
 Subject: reconstruct waveformm from harmonic of FFT From: Matt J Date: 28 May, 2012 12:00:07 Message: 5 of 6 "Andrew Sun" wrote in message ... > > > 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).
 Subject: reconstruct waveformm from harmonic of FFT From: Andrew Sun Date: 28 May, 2012 12:29:05 Message: 6 of 6 "Matt J" wrote in message ... > "Andrew Sun" wrote in message ... > > > > > 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.