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
|
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.
//
// Copyright INRIA
[w,ind]=sort(-w);w=-w;mag=mag(ind);
dnum=length(w);
pw=sqrt(w(dnum)*w(1));
if pw<1.d-8 then error('invalid frequency range');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('w should be positive')
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('not sampled high enough')
end
fresp=fresp(:);
function m=fliplr(m)
//utility fct.
[p,q]=size(m);
m=m(:,q:-1:1);
|
