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.

I have a time series that I converted to the frequency domain using fft in matlab. I want to reduce the amplitude at a given frequency range (remove a peak between 1.5 and 2Hz) but keep the phase constant.

If y=fft(x,N) is my spectrum in the frequency domain, the amplitude is given by: ampl=abs(y) and phase=angle(y); I have calculated my corrected amplitude (amplnew), which looks reasonable when I plot frequency versus amplnew.
But when I try to rebuild y from the new amplitude (amplnew) and original phase (phase), I get large high frequency spikes in the frequency range of interest when I recalculate the amplitude from amp2=abs(ynew). It seems like I should be able to recalculate y using either: ynew(f)=amplnew(f)*exp(j*phase(f)) or ynew(f)=amplnew(f)*(cos(phase(f))+j*sin(phase(f))) where f is the frequency at each point in the frequency range of interest. But amp2=abs(ynew) looks nothing like amplnew. Can someone tell me what I am doing wrong? I am only making changes in a small frequency range and want the original amplitude and phase to stay the same at all other frequencies. (I tried post images but as this is my first post I don't have enough reputation)

Update: I realized that the "j" in the above equations for ynew(f) should be an "i". This fixes the high amplitude but not the oscillations. When plotted in the frequency domain, the amplitude abs(ynew) should decrease linearly between 1.5 and 2 Hz. Instead I see a sinusoid in that range where the amplitude of the peaks of the sinusoids would match the amplitude I want (amplnew).

Update: Sorry for the formatting above. If it helps, here is the relevant parts of the actual code

 %Getting the time series data  
 [headers,sacdata] = sacreada(file);  

%This is the time vector  
time=linspace(headers.b,headers.e,headers.npts);  

 dt=.1355;                 %time step  
 fs=1/dt;                  %frequency step  

 npts=headers.npts;        %original number of points  
 N2=ceil(log2(npts));  
 Nfft=2^N2;                %new number of points  

 y=fft(sacdata,Nfft);      %Taking the fft of the data  

 fnyq=fs/2;                %Nyquist frequency  
 Nfre=Nfft/2+1;  

 %This is the frequency  
 freq=linspace(0,fnyq,Nfre);  % FFT symmetric about the Nyquist frequency  

 %Calculating the amplitude and phase  
 ampl=abs(y);                 %Amplitude of y  
 phase=angle(y);              %Phase of y  

 %Finding the indexes of the frequencies in the range of interest  

 m=1;  
 for j=1:length(freq)  

   if(freq(j)>=1.5 && freq(j)<=2)  
     fixfreq(m)=j;  
     m=m+1;  
   end  

 end  

 ynew=y;  

 FL=length(fixfreq);        %Number of points to be fixed  

 ampl3=ampl;  

 %The next section of code looks at each frequency needing to be changed  
 % Corrected amplitude(ampl3) uses the equation of a line between the amplitudes at    
 % 1.5 and 2.0Hz  

 for k=1:FL  

   aref1=ampl(fixfreq(1));                           %Amplitude at 1.5Hz  
   fref1=freq(fixfreq(1));                           %Frequency at 1.5Hz  
   aref2=ampl(fixfreq(FL));                          %Amplitude at 2.0Hz  
   fref2=freq(fixfreq(FL));                          %Frequency at 2.0Hz  

   ms=(aref2-aref1)/(fref2-fref1);                   %Slope of the line  
   ampnew(k)=ms*freq(fixfreq(k))+(aref1-ms*fref1);   %Equation for a line: y=mx+b  

   ampl3(fixfreq(k))=ampnew(k);                      %This is the corrected amplitude  

   %Recreating y from amplitude and phase  
   ynew(fixfreq(k))=ampnew(k)*(cos(phase(fixfreq(k)))+i*sin(phase(fixfreq(k))));  

 end  

  %This should be the same as ampl3 (but isn't)  
  ampl2=abs(ynew);  

  %Taking the ifft to get the new time series  
 xnew=ifft(ynew);  
share|improve this question
    
The variables j and i have the same default value in Matlab, so that is not the problem. The methodology you describe sounds correct, but, honestly, your post is difficult to read due to bad formatting. Perhaps you can reformat it and use the preformatted text option to show the actual Matlab code you are using. It sounds like a programming error, but I can't help until I can see what exactly it is your are doing. –  AnonSubmitter85 Aug 15 '13 at 22:07
    
For some reason the plot did change which I changed j to i. I thought they should be the same also. I tried to put the matlab code I used into my post. I guess I should have tried the preformatted option... In any case, see if the Matlab code makes sense and clarifies what I'm doing –  Kristin Phillips Aug 15 '13 at 23:40
    
You are overwriting the default value of j when you use it in the loop, so that's why it isn't equal to $\sqrt{-1}$. (Matlab actually recommends that you use 1j for $\sqrt{-1}$.) I also just ran your code and get a maximum difference between ampl2 and ampl3 of 3.5527e-15, so I don't see where your problem is. –  AnonSubmitter85 Aug 16 '13 at 2:25
    
I found the issue. I had switched from using j as my imaginary number to i. But realized that my code had a loop that uses i as an index. So I switched my indexes to ii which fixed the issue. –  Kristin Phillips Aug 16 '13 at 17:35

1 Answer 1

Perhaps an example of what I think it is you are trying to do will help:

N = 1000;
x = rand(N,1) .* exp( 1j * 2*pi*rand(N,1) );  % Arbitrary input signal.
X = fft(x);                                   % Input spectrum.
Xnew = X;                                     % Copy the input.
newinds = 330:350;                            % Samples to be modified.
Xnew( newinds ) = 0.8 * Xnew( newinds );      % Attenuate by some amount.
xnew = ifft( Xnew );                          % Back to the time domain.

% Now we just plot the results.
subplot( 2,1,1 );
plot( abs( x ) );
hold on;
plot( abs( xnew ), 'r' );
title( 'Time Amplitude' );
subplot( 2,1,2 );
plot( abs( X ) );
hold on
plot( abs( Xnew ), 'r' );
title( 'Spectral Amplitude' );
share|improve this answer

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.