apsd_dpsd.m

来自「功率谱估计」· M 代码 · 共 64 行

%
function drms = APSD_DPSD(f,a)
%
n=length(f);
%
for i=1:n
    a(i)=a(i)*(386.^2)/((2.*pi*f(i))^4);
end
%
MAX = 5000;
%
ra=0.;
vrms=0.;
iflag=0;
%
for i=1:MAX
    s(i)=0.;
end
%
nn=length(f)-1;
%
for  i=1:nn
    if(f(i) < .0001)
        f(i)=.0001;
    end
%
    if(  f(i) <=0 )
        disp(' frequency error ')
        out=sprintf(' f(%d) = %6.2f ',i,f(i));
        disp(out)
        iflag=1;
    end
    if(  a(i) <=0 )
        disp(' amplitude error ')
        out=sprintf(' a(%d) = %6.2f ',i,a(i));
        disp(out)
        iflag=1;
    end  
    if(  f(i+1) < f(i) )
        disp(' frequency error ')
        iflag=1;
    end  
    if(  iflag==1)
        break;
    end
%    
    s(i)=log10( a(i+1)/ a(i) )/log10( f(i+1)/f(i) );
%   
 end
 %
 %disp(' RMS calculation ');
 %
 if( iflag==0)
    for i=1:nn
        if(s(i) < -1.0001 |  s(i) > -0.9999 )
            ra = ra + ( a(i+1) * f(i+1)- a(i)*f(i))/( s(i)+1.);
        else
            ra = ra + a(i)*f(i)*log( f(i+1)/f(i));
        end
    end
    drms=sqrt(ra);
end