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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA - 1988 - C. Bunks
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
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 :
// - if vector, the first signal data
// - amount of input data if scalar. In this case the first signal
// is read by batch by a user defined function xbatch=getx(n,offset).)
// y : optional, if present cross correlation is made else auto correlation
// - if vector, the second signal data .It must have the same number of element than x
// - if scalar, the second signal
// is read by batch by a user defined function ybatch=gety(n,offset).)
// sec_step : offset of each data window.
// if sec_step==sec_leng/2 50% overlap is made
// sec_leng : Number of points of the 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
//!
[lhs,rhs]=argn(0);
if sec_step>=sec_leng then
error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Argument #%d expected to be less than argument #%d.\n"),"pspect",1,2,1,2));
end
if and(wtype<>["re","tr","hm","hn","kr","ch"]) then
error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
"pspect",3,"''re'',''tr'',''hm'',''hn'',''kr'',''ch''"));
end
//Analyze calling sequence and construct window
if rhs==4 then,
//pspect(sec_step,sec_leng,wtype,x)
w=window(wtype,sec_leng);
crossFlag=%f; //autocorrelation
elseif rhs==5 then,
//pspect(sec_step,sec_leng,wtype,x,y) or pspect(sec_step,sec_leng,wtype,x,wpar)
if wtype=="kr" then,//pspect(sec_step,sec_leng,'kr',x,wpar)
wpar=y;
w=window(wtype,sec_leng,wpar);cwp=[];
crossFlag=%f; //autocorrelation
elseif wtype=="ch" then,//pspect(sec_step,sec_leng,'ch',x,wpar)
wpar=y;
[w,cwp]=window(wtype,sec_leng,wpar);
crossFlag=%f; //autocorrelation
else,//pspect(sec_step,sec_leng,wtype,x,y)
w=window(wtype,sec_leng);cwp=[];
crossFlag=%t; //cross-correlation
end,
else //pspect(sec_step,sec_leng,wtype,x,y,wpar)
[w,cwp]=window(wtype,sec_leng,wpar);
crossFlag=%t; //cross-correlation
end,
wpower=w*w';//the window energy
//Make x and y row vectors
x=matrix(x,1,-1);
if crossFlag then
y=matrix(y,1,-1);
if size(x,"*")<>size(y,"*") then
error(msprintf(gettext("%s: Arguments #%d and #%d must have the same sizes.\n"),"pspect",4,5));
end
end
//Average periodograms
sm=0*w;
if size(x,"*")==1 then //Batch processing of x and y data
//get number of sections to be used
nsecs=int((x-sec_leng+sec_step)/sec_step);//x contains the number of points of the signals
ovrlp=sec_leng-sec_step;
xd=[0*ones(1,sec_step) getx(ovrlp,1)];
if crossFlag 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));
fx=fft(xw,-1);
if crossFlag 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));
sm=sm+fx.*conj(fft(yw,-1));
else
sm=sm+real(fx.*conj(fx));
end
end
else // Signal data given in x and y if cross-correlation
//get number of sections to be used
nsecs=int((size(x,"*")-sec_leng+sec_step)/sec_step);
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));
fx=fft(xe,-1);
if crossFlag then
yd=y(indd);
ye=w.*(yd-(sum(yd)/sec_leng));
sm=sm+fx.*conj(fft(ye,-1));
else
sm=sm+real(fx.*conj(fx));
end
end
end
//Normalization
sm=sm/(nsecs*wpower);
endfunction
|
