Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Rebuilding original signal from frequencies, amplitude, and phase obtained after doing an fft.

Greetings

I'm trying to rebuild a signal from the frequency, amplitude, and phase obtained after I do an fft of a signal, but when I try and combine the fft data (frequency, amplitude, and phase) back to see if I get a similar signal, the pattern is a little off. I think it has to do with my formula which may be a little incorrect.

The formula I'm using to combine the data is:

amplitude*sin(2*pi*time*frequency+phase)+amplitude*cos(2*pi*time*frequency+phase);

Please note: At the moment I don't want to use IFFT due to the fact that I will be editing the amplitude and frequencies before the calculations/rebuilding of the signal are done

An image of the plot is below. The top one is the original signal and the bottom one is the signal created with the equation. If you want to know I'm using matlab but I think the problem is with the equation. enter image description here


Here's the function I use to get the frequency, amplitude and phase of the signal.

function [freq,amp,phase,sigfft,sigifft,sigphase]=rtplotfft(vp_sig_orig,Fs)
    vp_sig_orig=vp_sig_orig';
    vp_sig_len=length(vp_sig_orig); %get sample rate from vp fs_rate needs to be an even number?

    % Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
    nfft= 2^(nextpow2(length(vp_sig_orig))); 
    % Take fft, padding with zeros so that length(fftx) is equal to nfft 
    fftx = fft(vp_sig_orig,nfft); 
    sigfft= fft(vp_sig_orig); 
    sigifft=ifft(sigfft);
    sigphase = unwrap(angle(sigfft')); %get phase of orginal signal
    % Calculate the number of unique points
    NumUniquePts = ceil((nfft+1)/2); 
    % FFT is symmetric, throw away second half 
    fftx = fftx(1:NumUniquePts); 
    % Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
    mx = abs(fftx)/length(vp_sig_orig); 

    % Take the square of the magnitude of fft of x. 
    %mx = mx.^2; rem'd out to get amplitude to work

    % Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
    % The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.
    if rem(nfft, 2) % odd nfft excludes Nyquist point
    mx(2:end) = mx(2:end)*2;
    else
    mx(2:end -1) = mx(2:end -1)*2;
    end

    %yamp=(mx(1,:)/max(abs(mx(1,:)))*1); % keep at 1, dont use in function use in script file directly
    amp=mx;

    % This is an evenly spaced frequency vector with NumUniquePts points. 
    freq_vect = (0:NumUniquePts-1)*vp_sig_len/nfft; 
    freq=freq_vect';

    %get phase of new signal
    phase = unwrap(angle(fftx)); %get phase of orginal signal 

tia

share|improve this question
    
Is there a question here? –  Doug Spoonwood Jun 2 '11 at 1:31
1  
Is there a reason you're using both a sine term and a cosine term? Since you're already keeping track of phase, I would expect you to use either sine or cosine, but not both. –  Jim Belk Jun 2 '11 at 3:11
    
You're processing the phase and amplitude differently. Is there a reason for that? The sigphase line is operating on the original-size fft, but the mx line is operating on the power-of-2 fft? You're throwing away the redundant second half of the magnitude information, but not of the phase information? You're getting the magnitude information from the power-of-2 fft but then dividing by the original-size length to normalize it? Wouldn't it be something like mx = abs(fftx)/(nfft/2);? –  endolith Jun 8 '11 at 19:03
    
What's the point of sigifft=ifft(sigfft)? sigifft should be identical to vp_sig_orig. Is it? –  endolith Jun 8 '11 at 19:07

1 Answer 1

up vote 1 down vote accepted

This is what the IFFT is for. Do your modifications to the spectrum, and then get them back into the form IFFT expects, and then do an IFFT. Don't bother trying to re-implement the IFFT when you don't have to. It will be slower and more error-prone.


So each point in an FFT is a complex number that contains information about the phase and amplitude at that frequency. You're extracting this information with these lines:

sigphase = unwrap(angle(sigfft')); %get phase of orginal signal
mx = abs(fftx)/length(vp_sig_orig); 

You're getting the magnitude and phase of each complex number, converting from rectangular complex number to polar form.

After modifying the phases and magnitudes, to reverse them back into something suitable for the IFFT, you just need to undo these and convert each back into a rectangular complex number, using the formula:

$z = r(\cos \varphi + i\sin \varphi ).\, $

where z is the resulting complex number, r is the magnitude, and φ is the phase. You'll have to undo the scaling of the magnitudes (and undo any change from radians to degrees, etc.) You'd do this for each mag, phase pair in the two vectors, combining them back into a single vector of the original length, made up of complex numbers. Then feed this to the IFFT. If you hadn't made any changes, it should return the original signal.

share|improve this answer
    
I can do the modifications but how do I get it back into a form IFFT can use? I have frequency, amplitude, and phase. –  Rick T Jun 8 '11 at 8:36
    
@Rick Do the opposite of what you did to convert them into frequency amplitude and phase. :) Can you detail your steps in the question? –  endolith Jun 8 '11 at 15:24
    
I'll include / add the function that creates the frequency, amplitude and phase to the question but I'm at a lost as how to reverse this process. –  Rick T Jun 8 '11 at 16:35
    
and everyone thanks again guys you all helped me so much it looks like I just needed to modify the equation a little. see code below: amplitude*cos(2*pi*time*frequency+phase) –  Rick T Jun 9 '11 at 17:09

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.