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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
|
function [sm,cwp]=pspect(sec_step,sec_leng,wtype,x,y,wpar)
//<sm,cwp>=pspect(sec_step,sec_leng,wtype,x,y,wpar)
//Cross-spectral estimate between x and y if both are given
//and auto-spectral estimate of x otherwise.
//Spectral estimate obtained using the modified periodogram method
// x :data if vector, amount of input data if scalar
// y :data if vector, amount of input data if scalar
// sec_step :offset of each data window
// sec_leng :length of each data window
// wtype :window type (re,tr,hm,hn,kr,ch)
// wpar :optional window parameters
// : for wtype='kr', wpar>0
// : for wtype='ch', 0<wpar(1)<.5, wpar(2)>0
// sm :power spectral estimate in the interval [0,1]
// cwp :unspecified Chebyshev window parameter
//!
//author: C. Bunks date: 14 Sept 1988
// Copyright INRIA
[lhs,rhs]=argn(0);
cross=0;
if sec_step<sec_leng then,
//get number of sections to be used
xsize=maxi(size(x));
if xsize==1 then,
xsize=x;
end,
nsecs=int((xsize-sec_leng)/sec_step);
//construct window
if rhs==4 then,
w=window(wtype,sec_leng);
else if rhs==5 then,
if wtype=='kr' then,
wpar=y;
w=window(wtype,sec_leng,wpar);
else if wtype=='ch' then,
wpar=y;
[w,cwp]=window(wtype,sec_leng,wpar);
else,
cross=1;
w=window(wtype,sec_leng);
end,
end,
else,
[w,cwp]=window(wtype,sec_leng,wpar);
cross=1;
end,
end,
wpower=w*w';
//average periodograms
sm=0*w;
if maxi(size(x))==1 then,
ovrlp=sec_leng-sec_step;
xd=[0*ones(1,sec_step) getx(ovrlp,1)];
if cross==1 then,
yd=[0*ones(1,sec_step) gety(ovrlp,1)];
end,
for k=1:nsecs,
xd(1:ovrlp)=xd(sec_step+1:sec_leng);
xd(ovrlp+1:sec_leng)=...
getx(sec_step,sec_leng+(k-1)*sec_step+1);
xw=w.*(xd-(sum(xd)/sec_leng)*ones(1:sec_leng));
fx=fft(xw,-1);
if cross==1 then,
yd(1:ovrlp)=yd(sec_step+1:sec_leng);
yd(ovrlp+1:sec_leng)=...
gety(sec_step,sec_leng+(k-1)*sec_step+1);
yw=w.*(yd-(sum(yd)/sec_leng)*ones(1:sec_leng));
fy=fft(yw,-1);
sm=sm+fx.*conj(fy);
else,
sm=sm+real(fx.*conj(fx));
end,
end,
else,
ind=(1:sec_leng);
for k=1:nsecs,
indd=ind+(k-1)*sec_step*ones(1:sec_leng);
xd=x(indd);
xe=w.*(xd-(sum(xd)/sec_leng)*ones(1:sec_leng));
fx=fft(xe,-1);
if cross==1 then,
yd=y(indd);
ye=w.*(yd-(sum(yd)/sec_leng)*ones(1:sec_leng));
fy=fft(ye,-1);
sm=sm+fx.*conj(fy);
else,
sm=sm+real(fx.*conj(fx));
end,
end,
end,
sm=sm/(nsecs*wpower);
//input error
else,
error('section step is bigger then section length');
end
|
