1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
|
////////////////////////
//SPECTRAL ESTIMATION///
////////////////////////
// Copyright INRIA
xselect();xbasc();
//generate white data
rand('normal');
rand('seed',0);
x=rand(1:1024-33+1);
//make low-pass filter with eqfir
nf=33;
bedge=[0 .1;.125 .5];
des=[1 0];
wate=[1 1];
h=eqfir(nf,bedge,des,wate);
//filter white data to obtain colored data
h1=[h 0*ones(1:maxi(size(x))-1)];
x1=[x 0*ones(1:maxi(size(h))-1)];
hf=fft(h1,-1);
xf=fft(x1,-1);
yf=hf.*xf;
y=real(fft(yf,1));
//plot magnitude of filter
h2=[h 0*ones(1:167)];
hf2=fft(h2,-1);
hf2=real(hf2.*conj(hf2));
hsize=maxi(size(hf2));
fr=(1:hsize)/hsize;
plot2d(fr',log(hf2)');
xtitle('Data spectrum','frequency','magnitude');
halt();xbasc();
//pspect example
[sm1]=pspect(100,200,'tr',y);
smsize=maxi(size(sm1));
fr=(1:smsize)/smsize;
plot2d(fr',log(sm1)')
xtitle('Spectral estimation','frequency','spectral power');
halt();xbasc();
//cspect example
[sm2]=cspect(100,200,'tr',y);
smsize=maxi(size(sm2));
fr=(1:smsize)/smsize;
plot2d(fr',log(sm2)');
xtitle(['Spectral estimation ; periodogram method'],' ' ,' ' )
halt();
xbasc();
|
