HW.rar_harmonic filter_小波滤波器_带通滤波_谐波小波_谐波滤波器

clc; clear all; fs=5000; N=10000; fid=fopen('c:\yibanEMD.txt');%%%%%读数据 x=fscanf(fid,'%f',N); f1=330; f2=360; % prepare the signal xsize=size(x); if xsize(2)==1 x=x'; end L=length(x); df=fs/L; if nargin<5 bf=f2-f1;% bandwidth if bf<1 df=min(bf/10,df); end N=ceil(fs/df);% Nfft end X=fft(x,N); % N points FFT of the signal want=abs(X);%看原信号频谱，便于自己理解 want1=ifft(X,N);%看逆傅里叶变换之后时域信号，便于自己理解 df1=floor(f1/df)+1; % convert the analogue frequency to digital form df2=floor(f2/df)+1; W=zeros(1,N); XW=zeros(1,N); for i=df1:df2 W(i)=1;% wavelet function in frequency domain, only the desired band nonzero end want2=conj(W);%看结果，便于自己理解 XW=X.*conj(W);% multiplication in frequency domain want3=abs(XW);%看HWT之后频谱，便于自己理解 x_hwt=ifft(XW,L); want4=x_hwt';%看结果，便于自己理解 %RDT N=length(x_hwt); n=9000; s=1.5*std(x_hwt); l=[]; for i=1:N-1 y=x_hwt-s; if sign(y(i))~=sign(y(i+1)) l=[l i]; end if i>=N-n; break; end end y=[]; for i=1:n u=0; for j=1:length(l) u=u+x_hwt(l(j)+i-1); end y=[y u/length(l)]; end frequency(y,fs);%频谱 axis tight; % % plot wavelet % figure,subplot(211),plot([0:df:fs/2-df],W([1:floor(fs/2/df)])), % xlabel('frequency(Hz)'),title('Harmonic wavelet'); % t=[-4+1/fs:1/fs:4]; % Lt=length(t); % w=ifft(W,Lt);% time signal % w2=[conj(w(4*fs:-1:1)),w(1:4*fs)]; % subplot(212),plot(t,real(w2),'r'),hold on,plot(t,imag(w2)), % xlabel('time(s)');