psd_icf.m

能在一些信号处理中用到的

function  pis = psd_icf(s,Fs,Bd,pt)
% s      input signal
% Fs     sample frequency
% Bd     带通滤波器带宽
% de_time      delay_time
% N = length(s);
% 列向量


de_time0 = fix(Fs/(Bd/2));
% dt = fix(de_time0/6);
% 
%  dt_sequ = [3 ]*dt;
dt_sequ = fix([0.35 0.4 0.45]*de_time0);
%  dt_sequ = [ 40  60 80];
% dt_sequ = fix(de_time0/2);


s = s(:);
% N = 2*pow2(nextpow2(length(s)));

Np = pow2(7);
L = Np/4;
df  = Fs/Np;
M = length(s)/Np;
dalpha = df/M;
P = pow2(nextpow2(Fs/dalpha/L));
N = P*L;


N2 = fix(N/2);
B = 2*fix(Bd*N/Fs/4);
kk = ((-N/2):(N/2-1))*Fs/N;
kk0 = wkeep(kk,N2 - 200,'l');
sfm = zeros(N,1); 


    Xr =s.*wshift('1d',s,0);
    b = fft(Xr,N);c = fftshift(b);
    c = c.*conj(c)/N;
    c = c/max(c);
    cf = wkeep(c,N2 - 200,'l');
    [zh_max,pisi] = max(abs(cf));
    pis_fc = abs(kk0(pisi))/2;
if(pt)
   figure;
    plot(kk, c.^(1/1));grid on;
    title('PSD delay time = 0');
    xlabel ('frequency'); ylabel('PSD magnitude');
end
   
if (pisi-2*B <= 0)
    pisi = pisi + 2*B +10;
end
 
% Xr =s.*wshift('1d',s,68);
%       b = fft(Xr,N);c = fftshift(b);
%       c = c.*conj(c)/N;
%       c = c/max(c);
%       figure;   plot(kk,c);grid on 
%       title(['PSD delay time = ',int2str(68)]);
%       xlabel ('frequency'); ylabel('PSD magnitude');


  for i = 1:length(dt_sequ)  
%        subplot(3,1,i);
      Xr =s.*wshift('1d',s,dt_sequ(i));
      b = fft(Xr,N);c = fftshift(b);
      c = c.*conj(c)/N;
      c = c/max(c);
      sfm = sfm + c;
%       figure;   plot(kk,c);grid on
%       for j = 1:30
%           [zh_max,pisi_max] = max(c);
%           c(pisi_max)=0;
%       end
 yc = wkeep(c,4*B,pisi-2*B);
 ykk = wkeep(kk,4*B,pisi-2*B);
 
%   yc0 = wkeep(c,8*B+1);
%  ykk0 = wkeep(kk,8*B+1);
%  figure;plot(ykk0,yc0);grid on
%  yc = wkeep(c,6*B,pisi-3*B);
%  ykk = wkeep(kk,6*B,pisi-3*B);

%  figure;plot(ykk,yc);grid on
%  dt = dt_sequ(i);
%       title(['PSD delay time = ',int2str(dt_sequ(i))]);
%       xlabel ('frequency'); ylabel('PSD magnitude around f = -2fc ');
% 
%       figure;   plot(kk,c);grid on 
%       title(['PSD delay time = ',int2str(dt_sequ(i))]);
%       xlabel ('frequency'); ylabel('PSD magnitude');
      
   end
%    subplot(4,1,4);
sfm = sfm/max(sfm);
% if(pt)
% %     figure;
% %      plot(kk,sfm.^3);grid on;
%      figure;    plot(kk,sfm.^1);grid on;
%      xlabel ('frequency'); ylabel('PSD magnitude');
% end
if pisi- 2*B <=0 
    pisi = 2*B + 10;
end
 ysfm = wkeep(sfm,4*B,pisi-2*B);
 ykk = wkeep(kk,4*B,pisi-2*B);
%   figure;plot(ykk,ysfm)
 
 k3 = [ones(1,1.5*B),zeros(1,B),ones(1,1.5*B)]';
 k3 = wkeep(k3,length(ysfm));
 ys3 = ysfm.*k3;
 [zh_max,pisi_max] = max(ys3);
 pis_fd = abs(ykk(pisi_max) - kk0(pisi));

 pis = [pis_fc; pis_fd];