pade.m

pade近似法频谱分析

function spec=pade(data,dt,f)
fre=f*1e9;
hsf=1/dt/2;
temp=0.8*hsf/max(fre);
decirate=fix(log(temp)/log(2));
for di=1:decirate
    data=decimate(data,2);
end
[siz1,siz2]=size(data);
if siz1==1
    N=siz2;
    if mod(N,2)==0
        N=N-1;
    end
    u=data(1:N).';
elseif siz2==1
    N=siz1;
    if mod(N,2)==0
        N=N-1;
    end
    u=data(1:N);
end

tv=(0:N-1).';
spec=zeros(size(fre));
num=max(size(fre));
for di=1:num
    cn=u.*exp(-i*2*pi*fre(di)*(2^decirate*dt)*tv);
    spec(di)=fastpade(cn);
end
spec=abs(spec.^2);
spec=spec/max(spec);
plot(f,spec);
xlabel('Frequency (GHz)');
ylabel('Intensity (a.u.)');
%---------------------founction fasepade-----------------------------------
function y=fastpade(c)
N=size(c,1);
cyita0=zeros(N,1);
cyita1=zeros(N,1);
czyita1=zeros(N,1);
cyita0=c;
cyita1=[c(1:N-1);0];
czyita1(2:N)=cyita1(1:N-1);
cita0=1;
cita1=1;
for j=1:(N-1)/2
    cita0=cita0-cyita0(N-j+1)/cyita1(N-j)*cita1;
    cyita0=cyita0-cyita0(N-j+1)/cyita1(N-j)*czyita1;
    cita1=(cyita0(N-j)*cita1-cyita1(N-j)*cita0)/(cyita0(N-j)-cyita1(N-j));
    cyita1=(cyita0(N-j)*cyita1-cyita1(N-j)*cyita0)/(cyita0(N-j)-cyita1(N-j));
    czyita1(2:N)=cyita1(1:N-1);
end
y=sum(cyita0)/cita0;