Prev: Sample time - sequence stair
Next: Help test tomorrow in the MATLAB
From: dbd on 7 Jun 2010 19:56
On Jun 7, 1:17 pm, "Erwin " <car...(a)gmail.com> wrote:
> dbd <d...(a)ieee.org> wrote in message <60c6c697-8145-4eb2-8f49-c89c2f06e....(a)z15g2000prh.googlegroups.com>...

> ,,,
> > When you have gotten this far, you can implement the last two dashed/
> > indented points in the first column on page 709 in your paper by:
>
> > 1) allocate a length 2*N+2 array of zeros
>
> > 2) fill the array with Y starting with the second element going
> > forward
>
> > 3) fill the array from the last element forward with the complex
> > conjugate of Y (this is where ifft expects the negative frequency data
> > to be, in reverse order)
>
> > 4)  ifft the array
>
> > 5) discard the imaginary component
> ...
> > Dale B. Dalrymple
>
> Hello Dale,
>
> Let me see if I understood what you've said.
>
> So I should create a new array, starting with a zero, followed by the values I created, and after them, I plug the complex conjugate of these values, and end with a zero?
>
> After I do the ifft, the size of the array is the (double + 2) of the one I created before, then how to I create the time function now? There are still imaginary components, but they aren't close to zero..
>
> Thanks a lot!

Note that I said "when you have gotten this far". Before you get there
you need to fix some things like your calculation of the array f. It
should have values greater than 0 and no greater than FS. There may be
other errors, I can't run your code because I don't have the functions
you call.

And to clean up my posting:
Make the new array size 2*N
What goes in the new array is a zero, Y(1:N-1) in order, a zero and
then conj(Y((N-1):-1:1) (note the reverse order), then ifft the array.
The leading zero is at DC and replaces the infinite value the
distribution would require there. Your simulated data is only valid
within a portion of the band from DC to Fs/2.

Dale B. Dalrymple
From: Erwin on 7 Jun 2010 21:04
Greg Heath <heath(a)alumni.brown.edu> wrote in message <87c4034a-9337-453e-a1e2-767ba19dc253(a)z8g2000yqz.googlegroups.com>...
> On Jun 7, 12:26 pm, "Erwin " <car...(a)gmail.com> wrote:
> > TideMan <mul...(a)gmail.com> wrote in message <b2eb37a4-0b03-4f38-a284-86738ae54...(a)z15g2000prh.googlegroups.com>...
> > > This is NOT the way to do it.
> > > Search on "phase randomisation" in this forum and you'll get an
> > > algorithm that does exactly what you want, but in the proper
> > > (numerically robust) way.
> >
> > Hello,
> >
> > Thanks for your reply!
> > Well, I've searched this forum and couldn't really find a response to what i want to implement...
> > I need to generate a Power Spectrum which follows a f^(-beta) law, and then have this in the time domain.
> >
> > What's really wrong with the code I've put in a previous message?
> >
> > Thanks again for your help, you're all very kind!
>
> The power law you state is just an asymptotic approximation
> that is not valid near zero frequency, Therefore, you
> should not just blindly use ifft and expect to find a
> discrete sampling of a nonpathological square integrable
> time function.
>
> Why not try something like (-pi*Fs <= w=2*pi*f < pi*Fs)
>
> abs(X) = a^(b/2)/(w^2+ a^2)^(b/4), b > 1.
>
> Although specified probability distributions could be
> generated, is is sufficient, without those specifications,
> to just choose a phase that is uniformly distributed
> over [0, 2*pi].
>
> Hope this helps,
>
> Greg
>

Hello Greg,

I didn't quite understood what you've said...
I'm going to post the code I got until now in the next response to Dale, could you give me your opinion?

Thanks a lot.
From: Erwin on 7 Jun 2010 21:49
> Note that I said "when you have gotten this far". Before you get there
> you need to fix some things like your calculation of the array f. It
> should have values greater than 0 and no greater than FS. There may be
> other errors, I can't run your code because I don't have the functions
> you call.
>
> And to clean up my posting:
> Make the new array size 2*N
> What goes in the new array is a zero, Y(1:N-1) in order, a zero and
> then conj(Y((N-1):-1:1) (note the reverse order), then ifft the array.
> The leading zero is at DC and replaces the infinite value the
> distribution would require there. Your simulated data is only valid
> within a portion of the band from DC to Fs/2.
>
> Dale B. Dalrymple

Hi Dale,

I want to generate a similar power law, which is between the low and high frequency I define below. The way you said, I get an entire array of real values after de ifft!
Here's the entire code:

clear all;
close all;

beta=0.5;
a=500;

Hf = 500; %high frequency
Lf = 0.03; %low frequency
Fs = 1/5; %sampling frequency
f = Lf:Fs:Hf; %frequency array
L = length(f);

for i = 1:L,
S(i) = (1/f(i))^(beta); %power function S(w)=(1/w)^beta
r(i) = a.*normrnd(0,(1/2)*S(i))*sqrt((1/2)*S(i)); %real part
im(i) = a.*normrnd(0,(1/2)*S(i))*sqrt((1/2)*S(i)); %imaginary part
G(i) = complex(r(i),im(i)); %complex function
end

Ytop = [0,G(1:L-1)]; %top part to plug into the ifft
Yconj = [0,conj(G(L-1:-1:1))]; %bottom part to plug into the ifft
Y = [Ytop,Yconj]; %merge of top and bottom parts

x = ifft(Y); %inv. FT
t = (0:2*L-1); %time array

Power = abs(Y).^2/(2*L); %power spectrum

areax = sum(abs(x).^2);
areaPower = sum(Power);
areaRatio=areax/areaPower; % parseval theorem verification

subplot(2,1,1);
loglog(f,Power(1:L)); %power spectrum plot
grid on;
grid minor;
title('|F(y)|^2')

subplot(2,1,2);plot(t,real(x)); %time signal plot
grid on;
grid minor;
title('f(x)');

Tell me what you think of the results.

Thank you again!
From: TideMan on 7 Jun 2010 22:25
On Jun 8, 1:49 pm, "Erwin " <car...(a)gmail.com> wrote:
> > Note that I said "when you have gotten this far". Before you get there
> > you need to fix some things like your calculation of the array f. It
> > should have values greater than 0 and no greater than FS. There may be
> > other errors, I can't run your code because I don't have the functions
> > you call.
>
> > And to clean up my posting:
> > Make the new array size 2*N
> > What goes in the new array is a zero, Y(1:N-1) in order, a zero and
> > then conj(Y((N-1):-1:1) (note the reverse order), then ifft the array.
> > The leading zero is at DC and replaces the infinite value the
> > distribution would require there. Your simulated data is only valid
> > within a portion of the band from DC to Fs/2.
>
> > Dale B. Dalrymple
>
> Hi Dale,
>
> I want to generate a similar power law, which is between the low and high frequency I define below. The way you said, I get an entire array of real values after de ifft!
> Here's the entire code:
>
> clear all;
> close all;
>
> beta=0.5;
> a=500;
>
> Hf = 500; %high frequency
> Lf = 0.03; %low frequency
> Fs = 1/5;  %sampling frequency
> f = Lf:Fs:Hf; %frequency array
> L = length(f);
>
> for i = 1:L,
>     S(i) = (1/f(i))^(beta); %power function S(w)=(1/w)^beta
>     r(i) = a.*normrnd(0,(1/2)*S(i))*sqrt((1/2)*S(i)); %real part
>     im(i) = a.*normrnd(0,(1/2)*S(i))*sqrt((1/2)*S(i)); %imaginary part
>     G(i) = complex(r(i),im(i)); %complex function
> end
>
> Ytop = [0,G(1:L-1)]; %top part to plug into the ifft
> Yconj = [0,conj(G(L-1:-1:1))]; %bottom part to plug into the ifft
> Y = [Ytop,Yconj]; %merge of top and bottom parts
>
> x = ifft(Y); %inv. FT
> t = (0:2*L-1); %time array
>
> Power = abs(Y).^2/(2*L); %power spectrum
>
> areax = sum(abs(x).^2);
> areaPower = sum(Power);
> areaRatio=areax/areaPower; % parseval theorem verification
>
> subplot(2,1,1);
> loglog(f,Power(1:L)); %power spectrum plot
> grid on;
> grid minor;
> title('|F(y)|^2')
>
> subplot(2,1,2);plot(t,real(x)); %time signal plot
> grid on;
> grid minor;
> title('f(x)');
>
> Tell me what you think of the results.
>
> Thank you again!

People have had trouble finding my phase randomisation thread (though
it comes up at the top of the list in Google Groups), but what it does
is to take the amplitude sqrt(S) and generate random phases between 0
and 2pi:
phase=rand(L,1)*2*pi; instead of real and imaginary parts using
randn.
Then, the Fourier transform is simply:
G=sqrt(S).*exp(i*phase);
etc

Note: you will have trouble here because you have used i as the index
in the for loop, which is a very silly thing to do when you are
dealing with complex numbers because it resets i from sqrt(-1) to some
other number.
Actually, though, you don't need a for loop at all.
Just use vectors.
From: Walter Roberson on 7 Jun 2010 22:33
TideMan wrote:
> generate random phases between 0
> and 2pi:
> phase=rand(L,1)*2*pi;

Unfortunately that will never generate a phase of 0 or 2 Pi -- though it
should get pretty close (2^(-23))
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4 5
Prev: Sample time - sequence stair
Next: Help test tomorrow in the MATLAB