Prev: Multiple Legends in GUI Axes
Next: resizing plots
From: martha on 14 Jul 2010 11:48
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
From: Wayne King on 14 Jul 2010 12:16
"martha " <dearmarylem(a)hotmail.com> wrote in message <i1km7k$bni$1(a)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 real-valued--you can't break the conjugate symmetry)
3.) Fill the vector of zeros with the complex-valued 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
From: martha on 14 Jul 2010 13:07
Thanx Wayne, thats an interesting point!
My aim is to get the mean absolute error between the initial signal and the estimator calculated with the most important frequencies. So i would actually need to approximate and compare the entire time series.
From: martha on 14 Jul 2010 13:39
I should maybe add that i have tried the code. it worked fine with the self defined series 'x' but unfortunately did not at all approximate the random time series i want to process.
From: Greg Heath on 14 Jul 2010 19:56
On Jul 14, 11:48 am, "martha " <dearmary...(a)hotmail.com> wrote:
> 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.

I am on a public machine that only contains IE and games. Therefore
I can only give you directions without details or valid code:

You can use the ifft:.

Y = fft(y);
PSD = abs(Y).^2/L; % Power Spectral density
P0 = sum(PSD) % proportional to the total power in y
P = cumsum(fliplr(sort(PSD))); % Cumulative sum of power
Find the index k that yields the % of power you wish to maintain
(say 90%)

P(k-1) < 0.9*P0 <= P(k)
Use the options in the sort command to find what frequencies {j}make
up
P(k).
Set all other values of Y equal to zero.
ifft the pruned Y.

Hope this helps.

Greg
 | 
Pages: 1
Prev: Multiple Legends in GUI Axes
Next: resizing plots