Phase correction in FFT

40 次查看（过去 30 天）
Jetson Ronald
Jetson Ronald 2013-5-12
Hello I have a signal (some harmonic function) and I want to increase the duration of signal using phase correction, i.e. adding phase to the existing phase. The reason why I use this approach is that these correction parameters are related to seismological parameters. In order to do that I followed the following steps.
1) Define the signal
2) Calculate FFT of the original signal, and take the amplitude and phase. 3) Define a frequency dependent function (envelope delay) for phase correction and integrate this function which is phase correction and add this correction to the existing Phase. I used a delta function.
4) Finally Build a new complex Fourier spectrum, using the corrected phase and the amplitude and take the inverse FFT to get the new signal.
The new signal should have increased duration, and the increased duration depends on the amplitude of the envelope delay function.
PROBLEM: My new signal has increased duration, I used 3 as amplitude of delta function and the duration of my signal is 5 sec and eventually my new signal has 8 sec but it has some more energy occurring after few seconds. Could someone help me fix the problem?.
Thank you
Jetson
Following is my code
if true
% %%%%Signal
Fs=50; % Sampling frequency
t_duration = 5; % Duration of signal
dt=1/Fs; % Sample time
t = 0:dt:t_duration; % Time vector
fsim=1; % Signal frequency
sn=sin(2*pi*fsim*t); % Signal
N=length(sn); % Length of signal
%%%% FFT of the Signal N1=round(2^(nextpow2(N)+1)); % No of frequencic, additioanl 1 is taken to account for the increased duration N12=N1/2; df=1/(dt*N1); % frequency interval Ny=N12*df; % Nyquist frequency f=0:df:Ny-df; % Frequency vector G=fft(sn,N1); % fft signal mag=abs(G); % amplitude of FFT of the signal phase=angle(G); % phase of the FFT of the signal
%%%%%%%%%%%%%%%%%%%%%%%%% % Phase Correction: FFT % %%%%%%%%%%%%%%%%%%%%%%%%%
%%% Delta Function for phase correction
fc=1; % Delta function has a value at this fequency amplitude_delta=3; % Amplitude of the delta function ff=0:df:fc; % dummy step for defining delta function funct1=zeros(1,length(ff)); % dummy step for defining delta function zz=length(f)-length(ff); % dummy step for defining delta function delta_fun=[funct1 amplitude_delta zeros(1,zz-1)]; % Delta Function (Envelope Delay) corec_ph=cumtrapz(f,(2*pi*delta_fun)); % Integral of Envelope Delay (Phase correction)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Amplitude Correction: FFT % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% amplitude of the modified signal is not changed so a unit amplitude % function is defined
amp_mag=1; corec_am=ones(1,length(f))*amp_mag; % Amplitude correction (no amplification)
%%%%%%%%%%%%%%%%% % FFT Operation % %%%%%%%%%%%%%%%%%
% initilizing Fourier componenets Amp=zeros(1,N1); % amplification for FFT manitude pha_cor=zeros(1,N1); % correction for Phase in FFT
for i=1:N12
Amp(i)=corec_am(i);
Amp(N1+1-i)=(Amp(i));
pha_cor(i)=-corec_ph(i);
pha_cor(N1+1-i)=-pha_cor(i);
end
%%%% Correction phase_cor=phase+pha_cor; % Modified Phase of the signal mag_amp=Amp.*mag; % Modified amplitude of the signal
%%%% Inverse FFT of Original signal R=mag.*cos(phase); I=mag.*sin(phase); spec=complex(R,I); sig_back=(ifft(spec,N1));
%%%% Corrected signal (Increased duration) R=mag_amp.*cos(phase_cor); I=mag_amp.*sin(phase_cor); spec=complex(R,I); sig_cor=real(ifft(spec,N1));
%%% Time vector of Corrected signal N=length(sig_back); T1=N*dt; Ts=[0:dt:(T1)-dt];
h1=figure(1); axes('position',[.1 .57 .38 .4]); stem(f,delta_fun,'r','linewidth',4,'markersize',2); hold on; set(gca,'fontsize',12) ylabel('Envelope Delay [sec]');grid on xlabel('Frequency [Hz]');axis tight
axes('position',[.57 .57 .38 .4]);
plot(f,-pha_cor(1:N12),'r','linewidth',4); hold on;
set(gca,'fontsize',12)
ylabel('Phase Correction [rad]');grid on
xlabel('Frequency [Hz]');
xlim([0 Ny]);
axes('position',[.1 .07 .85 .35]);
plot(Ts,sig_back,'b','linewidth',3); hold on; axis tight; grid on
plot(Ts,sig_cor,'r','linewidth',3); hold on; axis tight
set(gca,'fontsize',12);
ylabel('Amplitude');xlabel('Time [sec]');
legend('Original','Corrected'); legend boxoff
end
  1 个评论
Jetson Ronald
Jetson Ronald 2013-5-12
Sorry, the code I pasted above is not clear, so kindly refer this one instead of the above one.
%%%% Signal
Fs=50; % Sampling frequency
t_duration = 5; % Duration of signal
dt=1/Fs; % Sample time
t = 0:dt:t_duration; % Time vector
fsim=1; % Signal frequency
sn=sin(2*pi*fsim*t); % Signal
N=length(sn); % Length of signal
%%%%FFT of the Signal
N1=round(2^(nextpow2(N)+1)); % No of frequencic, additioanl 1 is taken to account for the increased duration
N12=N1/2;
df=1/(dt*N1); % frequency interval
Ny=N12*df; % Nyquist frequency
f=0:df:Ny-df; % Frequency vector
G=fft(sn,N1); % fft signal
mag=abs(G); % amplitude of FFT of the signal
phase=angle(G); % phase of the FFT of the signal
%%%%%%%%%%%%%%%%%%%%%%%%%
% Phase Correction: FFT %
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Delta Function for phase correction
fc=1; % Delta function has a value at this fequency
amplitude_delta=3; % Amplitude of the delta function
ff=0:df:fc; % dummy step for defining delta function
funct1=zeros(1,length(ff)); % dummy step for defining delta function
zz=length(f)-length(ff); % dummy step for defining delta function
delta_fun=[funct1 amplitude_delta zeros(1,zz-1)]; % Delta Function (Envelope Delay)
corec_ph=cumtrapz(f,(2*pi*delta_fun)); % Integral of Envelope Delay (Phase correction)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Amplitude Correction: FFT %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% amplitude of the modified signal is not changed so a unit amplitude
% function is defined
amp_mag=1;
corec_am=ones(1,length(f))*amp_mag; % Amplitude correction (no amplification)
%%%%%%%%%%%%%%%%%
% FFT Operation %
%%%%%%%%%%%%%%%%%
% initilizing Fourier componenets
Amp=zeros(1,N1); % amplification for FFT manitude
pha_cor=zeros(1,N1); % correction for Phase in FFT
for i=1:N12
Amp(i)=corec_am(i);
Amp(N1+1-i)=(Amp(i));
pha_cor(i)=-corec_ph(i);
pha_cor(N1+1-i)=-pha_cor(i);
end
%%%%Correction
phase_cor=phase+pha_cor; % Modified Phase of the signal
mag_amp=Amp.*mag; % Modified amplitude of the signal
%%%%Inverse FFT of Original signal
R=mag.*cos(phase);
I=mag.*sin(phase);
spec=complex(R,I);
sig_back=(ifft(spec,N1));
%%%%Corrected signal (Increased duration)
R=mag_amp.*cos(phase_cor);
I=mag_amp.*sin(phase_cor);
spec=complex(R,I);
sig_cor=real(ifft(spec,N1));
%%%Time vector of Corrected signal
N=length(sig_back);
T1=N*dt;
Ts=[0:dt:(T1)-dt];
h1=figure(1);
axes('position',[.1 .57 .38 .4]);
stem(f,delta_fun,'r','linewidth',4,'markersize',2); hold on;
set(gca,'fontsize',12)
ylabel('Envelope Delay [sec]');grid on
xlabel('Frequency [Hz]');axis tight
axes('position',[.57 .57 .38 .4]);
plot(f,-pha_cor(1:N12),'r','linewidth',4); hold on;
set(gca,'fontsize',12)
ylabel('Phase Correction [rad]');grid on
xlabel('Frequency [Hz]');
xlim([0 Ny]);
axes('position',[.1 .07 .85 .35]);
plot(Ts,sig_back,'b','linewidth',3); hold on; axis tight; grid on
plot(Ts,sig_cor,'r','linewidth',3); hold on; axis tight
set(gca,'fontsize',12);
ylabel('Amplitude');xlabel('Time [sec]');
legend('Original','Corrected'); legend boxoff

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Youssef Khmou
Youssef Khmou 2013-5-12
hi,
i suggest that this is due to Gibbs effect, i tried to change the amplitude, and time duration but always there is a small fluctuation around zero,
try :
t_duration = 20; % Duration of signal instead of 5 and l
  1 个评论
Jetson Ronald
Jetson Ronald 2013-5-13
Hello Youssef
Thank You for your answer, I am not sure this is a Gibbs effect, why do you think so? Do you see any other issue with the code? I tried it with duration 20 sec as you suggested still, I encounter the same problem.
I would very much apprciate your thoughts on this.
Thank you

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by