You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
106 lines
2.9 KiB
106 lines
2.9 KiB
function [sm,cwp]=specgram(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
|
|
endfunction
|
|
|