Phase correction in FFT

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

[Jetson Ronald]

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

Answer 1

[Youssef Khmou]

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

[Jetson Ronald]

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