dtfilter.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 先带通滤波，再小波包去噪
 function b = DTfilter(s,Fs,Fc,Bd)
% Fc 估计载频，Bd 带通滤波器带宽

if Fc >= Fs/2,b = s; return; end

s = s(:);
Nx = length(s);
% N = pow2(nextpow2(Nx));
N = Nx;
N2 = fix(N/2);   
B = fix(Bd*N/Fs/2);
kk = ((-N/2):(N/2-1))*Fs/N;
[zhi,pisiL]=min(abs(kk + Fc));

b = fft(s,N);c = fftshift(b);
c2 = c;
kk0 = 1:length(kk);

% if pisiL - B <= 0 | pisiL - B <= 0,b = s; return; end
% kk_fL  = wkeep(kk0,2*B,pisiL-B);c2(kk_fL) = 0;
% kk_fR  = wkeep(kk0,2*B,pisiR-B);c2(kk_fR) = 0;
if pisiL - B <= 0 
    pisiL = pisiL+B +10;
end

kk_fL = [pisiL-B:pisiL+B];  c2(kk_fL) = 0;
kk_fR = length(kk)-kk_fL+2;   c2(kk_fR) = 0;
c3 = c - c2;
bnew = ifftshift(c3);
b2 = ifft(bnew);
b = real(b2);

b = b(1:Nx);
%  a = fftshift(fft(b));
%  plot(abs(a))
%  a2 = conv(a,a);
%  plot(abs(a2))
% kk_new = (-N+1:N-1)*Fs/N; 
% kk_new = kk_new';
% kk2 = wkeep(kk_new,N+1);
% c4 = abs(c3);
% c5 = wkeep(c4,N+1);
% %  figure; plot(kk,c4);grid on
%  a2 = conv(c3,c3);
%  a4 = abs(a2);
%  a5 = wkeep(a4,N+1);
%  a5 = a5.^2;
%  figure;plot(kk2,a5/max(a5));grid on
%  figure;plot(real(c3))
%  figure;plot(imag(c3))
%  figure;plot(a4/max(a4))