"martha " <dearmarylem@hotmail.com> wrote in message <i1km7k$bni$1@fred.mathworks.com>...
> Hi all,
>
> i want to calculate the ifft of a transformed discrete time series. Since i only want to use the most important (largest) frequencies to calculate the inverse, I cannot use the function ifft. in my approximation code is, however, a magnitude shift  and i just cannot figure out how to fix it.
>
> Thanx for ur help!
> Martha
>
> % I have calculated the fft of a discrete sample set 'y' with length 'L'
>
> NFFT = 2^nextpow2(L);
> Y = fft(y,NFFT)/L;
> f = Fs/2*linspace(0,1,NFFT/2);
>
> LS = 2*abs(Y(1:NFFT/2));
>
> % Phase of FFT values
> for j=1:NFFT/2,
> Phi(j)=unwrap(angle(Y(j)));
> end
>
> % Sort according to importance of frequencies:
> [LS_sorted, idx] = sort(LS,'descend');
>
> % 'm' most important frequencies
> max_m = 64;
> MAE = NaN(1,max_m);
> for m=1:max_m,
>
> f_main = f(idx(1:m));
> LS_main = LS(idx(1:m));
> Phi_main = Phi(idx(1:m));
>
> data_approx = 0;
>
> for j=1:m,
> data_approx = data_approx + LS_main(j)*cos(2*pi*f_main(j)*t+Phi_main(j));
> % data_approx = data_approx + LS_main(j)*cos(Phi_main(j))*cos(2*pi*f_main(j)*t)  LS_main(j)*sin(Phi_main(j))*sin(2*pi*f_main(j)*t);
> end
Hi Martha, Out of curiosity why don't you just:
1.) Create a vector of zeros
2.) Find the largest magnitudes (make sure you take them in pairs if the signal is realvaluedyou can't break the conjugate symmetry)
3.) Fill the vector of zeros with the complexvalued coefficients you identified in step 2 (in the appropriate places of course!)
4.) Take the inverse Fourier transform of that.
For example:
Fs = 1000;
t = 0:(1/Fs):1(1/Fs);
x = cos(2*pi*100*t)+0.5*sin(2*pi*200*t)+0.25*randn(size(t));
xDFT = fft(x);
[freqsort, idx] = sort(abs(xDFT),'descend');
xrecon = zeros(size(xDFT));
% Take the 4 largest magnitudes you have two sinusoids but that makes for
% 4 complex exponentials
xrecon(idx(1:4))=xDFT(idx(1:4));
xrecon = ifft(xrecon,'symmetric');
% plot and compare focusing on the 1st 250 points
plot(xrecon(1:250));
hold on;
plot(x(1:250),'k');
This seems like what you are trying to do with your code anyway.
Hope that helps,
Wayne
