[Math] reducing amplitude of fft spectrum with constant phase

fourier analysisMATLAB

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);

Best Answer

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' );
Related Question