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
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA
//
// 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 fresp = cepstrum(w,mag)
// Uses the complex-cepstrum (Oppenheim & Schafer, Digital
//Signal Processing, p. 501) to generate, at the frequencies
//w, a complex frequency response fresp whose magnitude is
//equal to magnitude data mag and whose phase corresponds
//to a stable, minimum phase transfer function.
if ~isreal(w) then
error(msprintf(gettext("%s: Input argument #%d must be real.\n"),"cepstrum",1));
end
if ~isreal(mag) then
error(msprintf(gettext("%s: Input argument #%d must be real.\n"),"cepstrum",2));
end
[w,ind]=gsort(-w);w=-w;mag=mag(ind);
dnum=length(w);
pw=sqrt(w(dnum)*w(1));
if w(1)<=0 then
error(msprintf(gettext("%s: Wrong values for input argument #%d: Elements must be positive.\n"),"cepstrum",1));
end
if pw<1.d-8 then
error(msprintf(gettext("%s: Wrong values for input argument #%d: Distance between Max and Min elements must be greater than %g.\n"),..
"cepstrum",1,1d-8));
end
if size(w,"*")<>size(mag,"*") then
error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same numbers of elements expected.\n"),"cepstrum",1,2));
end
tcomes=(w.*w)/pw/pw;
dw=acos((1-tcomes)./(1+tcomes));
lowf=dw(1);
highf=dw(length(dw));
if lowf <= 0 | highf >=%pi
error(msprintf(gettext("%s: Wrong values for input argument #%d: Elements must be positive.\n"),"cepstrum",1));
end
nn=ceil((log(2*%pi/min(lowf,%pi-highf)))/log(2));
npts=2*(2^nn);hnpts =2^nn;
if npts<4096 then npts=4096;hnpts=2048;end
lmagin =(1/log(10))*log(mag);lindw =(%pi/hnpts)*(0:hnpts-1);
lmag=zeros(1,hnpts);topval=length(dw);
p=find(lindw<dw(1));lmag(p)=lmagin(1)*ones(1,length(p));
p=find(dw(topval)<=lindw);
lmag(p)=lmagin(topval)*ones(1,length(p));
for i=2:topval
p=find(lindw>=dw(i-1) & lindw<dw(i));
wrat=lindw(p) - dw(i-1)*ones(1,length(p));
wrat=(1/(dw(i)-dw(i-1)))*wrat;
lmag(p)=(ones(1,length(p))-wrat)*lmagin(i-1)+wrat*lmagin(i);
end
linmag=exp(log(10)*lmag);
dome=[lindw,(2*%pi)-fliplr(lindw)];
mag=[linmag,fliplr(linmag)];
ymag=log(mag.*mag);ycc=fft(ymag,1);nptso2=npts/2;xcc=ycc(1:nptso2);
xcc(1)=xcc(1)/2;xhat=exp(fft(xcc,-1));domeg=dome(1:2:nptso2-1);
xhat = xhat(1:nptso2/2);nptsslo=length(dw);nptsfst = length(domeg);
if domeg(1)<=dw(1) & domeg(nptsfst)>=dw(nptsslo)
fresp=zeros(1,nptsslo);
for i=1:nptsslo
p=min(find(domeg>=dw(i)));
wrat=(dw(i)-domeg(p-1))/(domeg(p)-domeg(p-1));
fresp(i)=wrat*xhat(p) + (1-wrat)*xhat(p-1);
end
else
error(msprintf(gettext("%s: Wrong values for input argument #%d: Not sampled high enough.\n"),"cepstrum",1));
end
fresp=fresp(:);
endfunction
function m=fliplr(m)
//utility fct.
[p,q]=size(m);
m=m(:,q:-1:1);
endfunction
|
