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