Prev: Sample time - sequence stair
Next: Help test tomorrow in the MATLAB
From: dbd on 7 Jun 2010 23:20
On Jun 7, 7:25 pm, TideMan <mul...(a)gmail.com> wrote:

> ...
> People have had trouble finding my phase randomisation thread
> ...

People have trouble finding your thread because you fail to provide an
explicit reference and because they fail to consistently misspell
randomization in searches.

Dale B. Dalrymple
From: Erwin on 7 Jun 2010 23:24
> 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.

Hello TideMan,

Yes, you're right I don't need the loop... it's just that I started with it. So I've followed what you've said, but when I plot the Power Spectrum I only get a line... why?

Here's what I've done:
clear all;

beta=2; %slope of loglog power
a=50; %amplitude

Hf = 10; %high frequency
Lf = 0.01; %low frequency
Fs = 0.001; %sampling frequency
f = Lf:Fs:Hf; %frequency array
L = length(f);

S = (1./f).^(beta);

phase=rand(1,L)*2*pi;
G=a.*sqrt(S).*exp(i.*phase);

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)');
From: dbd on 8 Jun 2010 00:28
On Jun 7, 8:24 pm, "Erwin " <car...(a)gmail.com> wrote:
> > 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.
>
> Hello TideMan,
>
> Yes, you're right I don't need the loop... it's just that I started with it. So I've followed what you've said, but when I plot the Power Spectrum I only get a line... why?
>
> Here's what I've done:
> clear all;
>
> beta=2; %slope of loglog power
> a=50; %amplitude
>
> Hf = 10; %high frequency
> Lf = 0.01; %low frequency
> Fs = 0.001;  %sampling frequency
> f = Lf:Fs:Hf; %frequency array
> L = length(f);
>
> S = (1./f).^(beta);
>
> phase=rand(1,L)*2*pi;
> G=a.*sqrt(S).*exp(i.*phase);
>
> 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)');

Erwin

I'm not going to waste time trying to keep up with a constantly moving
target when you fail to even clean up the problems people help you
with and persist in using function calls I don't have the toolbox for.
Here is an example of what I have said with vectors usage added and
Fortran-like indexing removed:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
beta=2; %slope of the log-log power law
a=500; %amplitude
N = 500; % half the number of samples
Fs = 1/5;
f = 2*(1:N)/Fs; % generated frequency space, includes no DC term

S = (1./f).^(beta); %power function S(w)=(1/w)^beta
Y = randn(1,N).*S + i*randn(1,N).*S;
newY = [0, Y(1:N-1) 0 conj(Y(N-1:-1:1))]; % As I said

figure
subplot(2,1,1)
loglog(abs(Y(1:N-1)))
grid on
title('Generated Spectrum Y')

timecx = ifft(newY);
sum(abs(real(timecx)))
sum(abs(imag(timecx)))
timerl = real(timecx);
ftest = abs(fft(timerl));

subplot(2,1,2)
loglog(abs(ftest(2:N-1)))
grid on
title('Spectrum calculated from time domain real samples')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

That leaves you with the exercises of determining unit scaling and of
controlling the bandwidth you wish to generate. Then you can try to
figure what you want about Parseval.

Dale B. Dalrymple
From: TideMan on 8 Jun 2010 00:45
On Jun 8, 3:20 pm, dbd <d...(a)ieee.org> wrote:
> On Jun 7, 7:25 pm, TideMan <mul...(a)gmail.com> wrote:
>
> > ...
> > People have had trouble finding my phase randomisation thread
> > ...
>
> People have trouble finding your thread because you fail to provide an
> explicit reference and because they fail to consistently misspell
> randomization in searches.
>
> Dale B. Dalrymple

Yes, it's an Antipodean practice that you need a licence to get the
spelling right when you are practising randomisation, especially with
coloured noise.
From: dbd on 8 Jun 2010 02:32
On Jun 7, 9:45 pm, TideMan <mul...(a)gmail.com> wrote:
> On Jun 8, 3:20 pm, dbd <d...(a)ieee.org> wrote:
>
> > On Jun 7, 7:25 pm, TideMan <mul...(a)gmail.com> wrote:
>
> > > ...
> > > People have had trouble finding my phase randomisation thread
> > > ...
>
> > People have trouble finding your thread because you fail to provide an
> > explicit reference and because they fail to consistently misspell
> > randomization in searches.
>
> > Dale B. Dalrymple
>
> Yes, it's an Antipodean practice that you need a licence to get the
> spelling right when you are practising randomisation, especially with
> coloured noise.

Yes, the US and the former British empire, two groups separated by
their common language.

Dale B. Dalrymple
First  |  Prev  |  Next  |  Last
Pages: 1 2 3 4 5
Prev: Sample time - sequence stair
Next: Help test tomorrow in the MATLAB