MATLAB Central Newsreader - Ifft on modified array fails

Ifft on modified array fails

Travis Bland — Fri, 06 Nov 2009 19:36:04 -0500

My end goal is to take an array, fft, set a range of frequencies equal to 0, and then ifft.
I get stuck once i set anything equal to 0. My ifft output is still complex. I stopped using my code, and just got an array in the workspace, manually changed a few entries to 0, and then tried plot(ifft(vf)) - i got scribbles, and complex results.

So i looked at some old code a friend wrote, he puts this in before he sets to 0. It's c++

for (int i = 0; i < nx; i++)
{
for (int j = 0; j < nyp; j++)
{
real(output1(i,j)) = real(output1(i,j)) /* sqrt(nx*ny)*/;
imag(output1(i,j)) = imag(output1(i,j)) /* sqrt(nx*ny)*/;
}
}

So i tried real(vf) = real(vf).... matlab didn't like that.

Any ideas?
Thanks,
Travis

Re: Ifft on modified array fails

Matt — Fri, 06 Nov 2009 21:10:21 -0500

"Travis Bland" <travisbland88@yahoo.com> wrote in message <hd1tr4$9lk$1@fred.mathworks.com>...
> My end goal is to take an array, fft, set a range of frequencies equal to 0, and then ifft.
> I get stuck once i set anything equal to 0. My ifft output is still complex. I stopped using my code, and just got an array in the workspace, manually changed a few entries to 0, and then tried plot(ifft(vf)) - i got scribbles, and complex results.
=======

We need to see more.

Why do you expect to get real output from the ifft? Is it a symmetric spectrum? Are the frequencies being set to 0 in such a way as to preserve this symmetry?

Re: Ifft on modified array fails

Travis Bland — Sat, 07 Nov 2009 04:01:02 -0500

"Matt " <xys@whatever.com> wrote in message <hd23bt$l4f$1@fred.mathworks.com>...
> "Travis Bland" <travisbland88@yahoo.com> wrote in message <hd1tr4$9lk$1@fred.mathworks.com>...
> > My end goal is to take an array, fft, set a range of frequencies equal to 0, and then ifft.
> > I get stuck once i set anything equal to 0. My ifft output is still complex. I stopped using my code, and just got an array in the workspace, manually changed a few entries to 0, and then tried plot(ifft(vf)) - i got scribbles, and complex results.
> =======
>
> We need to see more.
>
> Why do you expect to get real output from the ifft? Is it a symmetric spectrum? Are the frequencies being set to 0 in such a way as to preserve this symmetry?

Yes and yes. It is a symmetric spectrum and i believe i modified it to preserve that. It's data taken from ice samples in antarctica. And i guess i just assumed a complex answer was wrong, if i can just do abs(vt) and get good results, that's fine with me.

Here's some of my code, i just copied the parts for 0-100 mhz (and left 100-1000 out)

load std000.txt
vt=std000(:,2)

%power of vt over entire length of data file
for n = 1:length(vt) ;
vt_s = vt_s + (vt(n) .^2);
end

vf = (fft(vt));
vf_s =0;

%power of vf over entire file
for k= 1:length(vt);
    vf_s = vf_s + (abs( (vf(k)/sqrt(N)) .^2));
end

% 0 - 100 mhz removed
freq100 = vf;
bandwidth = 2000000000;
N = length(vt);
bin_ratio = bandwidth / N;
low_cut = 1 ; % user defined--- low cutoff
high_cut = 100 ; % user defined--- high cutoff
low_bin = (low_cut * 1000000) / ( bin_ratio);
high_bin = (high_cut * 1e6) / ( bin_ratio);

high_bin = ceil(high_bin);
low_bin = floor(low_bin);

for i = 1:high_bin;
    freq100(i) = 0;
end

for p = (N - high_bin):(N - low_bin);
    freq100(p) = 0;
end

%skipping through same calculations for different frequency ranges... we get to the %summation section

vt_100 = ifft(freq100);
vt_100_s = 0;
for n = 1:length(vt) ;
vt_100_s = vt_100_s + (vt_100(n) .^2);
end

Thanks!

Re: Ifft on modified array fails

Matt — Sat, 07 Nov 2009 15:12:01 -0500

"Travis Bland" <travisbland88@yahoo.com> wrote in message <hd2rdu$l9g$1@fred.mathworks.com>...

> high_bin = ceil(high_bin);
> low_bin = floor(low_bin);
>
>
>
> for i = 1:high_bin;
> freq100(i) = 0;
> end
>
>
> for p = (N - high_bin):(N - low_bin);
> freq100(p) = 0;
> end
>

I don't see anything symmetric in this truncation.

The important thing to remember is that, for a given index k, if you set
freq100(k)=0, then you must also set
freq100(N+1-k)=0 to preserve conjugate symmetry of the spectrum.

So, if you set freq100(1:high_bin)=0, as above, you must also set
freq100(N+1-(1:high_bin))=0 as well

Re: Ifft on modified array fails

Matt — Sat, 07 Nov 2009 16:04:01 -0500

"Matt " <xys@whatever.com> wrote in message <hd42o1$re2$1@fred.mathworks.com>...

> The important thing to remember is that, for a given index k, if you set
> freq100(k)=0, then you must also set
> freq100(N+1-k)=0 to preserve conjugate symmetry of the spectrum.

Make that freq100(N+2-k)=0

Re: Ifft on modified array fails

Travis Bland — Sun, 08 Nov 2009 08:17:03 -0500

"Matt " <xys@whatever.com> wrote in message <hd45ph$3ij$1@fred.mathworks.com>...
> "Matt " <xys@whatever.com> wrote in message <hd42o1$re2$1@fred.mathworks.com>...
>
> > The important thing to remember is that, for a given index k, if you set
> > freq100(k)=0, then you must also set
> > freq100(N+1-k)=0 to preserve conjugate symmetry of the spectrum.
>
> Make that freq100(N+2-k)=0

Well i could be mistaken but i think that's what i'm doing, i just do it in the second for loop. The plot is reflected at 1000mhz, when i run the code and remove 0-100mhz, 1900-2000mhz is also removed... so it stays symmetric, visually at least.

Is it deeper that? If i'm just wrong, i can accept that, take your advice and move on; i'm good at that! :)

thanks for your time!

Re: Ifft on modified array fails

Matt — Sun, 08 Nov 2009 16:27:01 -0500

"Travis Bland" <travisbland88@yahoo.com> wrote in message <hd5upv$nr4$1@fred.mathworks.com>...
> "Matt " <xys@whatever.com> wrote in message <hd45ph$3ij$1@fred.mathworks.com>...
> > "Matt " <xys@whatever.com> wrote in message <hd42o1$re2$1@fred.mathworks.com>...
> >
> > > The important thing to remember is that, for a given index k, if you set
> > > freq100(k)=0, then you must also set
> > > freq100(N+1-k)=0 to preserve conjugate symmetry of the spectrum.
> >
> > Make that freq100(N+2-k)=0
>
> Well i could be mistaken but i think that's what i'm doing, i just do it in the second for loop. The plot is reflected at 1000mhz, when i run the code and remove 0-100mhz, 1900-2000mhz is also removed... so it stays symmetric, visually at least.
===============

Well, visual tests can be misleading when the discrepancies are too small to show up noticably in a plot.

One clue here that something is wrong is that the number of steps in your first for-loop is high_bins, while the 2nd for-loop runs over
high_bins-low_bins+1 elements. So, the two loops are just not going to even be of the same length (let alone be symmetric) unless low_bins=1. If low_bins is equal to 1 all the time, I don't know why you would have it in there as a variable.

It occurs to me that maybe you want your first for-loop to start from low_bins instead of from 1, e.g.

for i = low_bin:high_bin;
freq100(i) = 0;
end

If so, you still need to make sure the 2nd for-loop uses the N+2-k rule, so it should be

range=low_bin:high_bin;

range(range==1)=[]; %No need to include DC

symrange=N+2-range; %The symmetric counterpart using the N+2-k rule

for p = symrange
freq100(p) = 0;
end

Re: Ifft on modified array fails

Matt — Sun, 08 Nov 2009 20:30:05 -0500

"Matt " <xys@whatever.com> wrote in message <hd6rgl$kgj$1@fred.mathworks.com>...

> for i = low_bin:high_bin;
> freq100(i) = 0;
> end
>
> If so, you still need to make sure the 2nd for-loop uses the N+2-k rule, so it should be
>
> range=low_bin:high_bin;
>
> range(range==1)=[]; %No need to include DC
>
> symrange=N+2-range; %The symmetric counterpart using the N+2-k rule
>
> for p = symrange
> freq100(p) = 0;
> end

Also, you would really do well to get acquainted with MATLAB's vector indexing operations and get rid of these for loops:

%For-loop Free !!!!

range=low_bin:high_bin;
symrange=N+2-range(range~=1);

freq100(range)=0;
freq100(symrange)=0;