齿轮振动信号的数字滤波处理-Matlab源代码

    clear all
    close all
    clc
    
    %% 读取齿轮振动信号
    x= load('齿轮振动信号.txt');     
    Fs=10240;      % 采样频率
    
    % 去除均值
    x=x-mean(x);
    xlen=max(size(x));
    if mod(xlen, 2)==1
    x=x(1:x1en-1);
    xlen= xlen-1
    end
    
    % handle1=figure;
    t=(0:xlen-1)/Fs;  %时间
    t=t';
    
    [a b]=size(x);
    if a==1
    x=x';
    end
    figure,plot(t, x);title('振动信号的波形图');
    xlabel('时间/s');
    ylabel('幅值/mv');
    grid on;
    % disp('按任意键继续…')
    % pause
    %hidem(handlel);
    
    %% 数字滤波处理
    
    disp('请用鼠标在谱图上选择滤波的截止频率');
    disp('左键选择,右键退出,按任意键继续')
    disp('1.低通滤波');
    disp('2.高通滤波')
    disp('3.带通滤波');
    filetype= input('请选择滤波的类型:');  % 输入1 表示低通；2表示选择高通 3 表示带通
    
    
    
    y=fft(x)/(xlen/2);     % 傅里叶变换
    absy=abs(y);
    
    angley=angle(y);
    z= abs(y(1:xlen/2));
    Amp=max(z);
    u=ones(1,xlen/2)* Amp;
    winlen= round(log2(xlen/40.));
    winlen=pow2(winlen);
    v=Amp*cos((-pi:(pi/winlen):0))/2+ Amp/2;
    w=Amp* cos((0:(pi/winlen):pi))/2+Amp/2;
    winlen= winlen+1;
    f=(0:(xlen-1)/2)*(Fs/xlen);
    Df=f(2);
    
    handle2=figure
    plot(f,z);
    title('信号的幅值谱图');
    xlabel('频率/Hz');
    ylabel('幅值/mv');
    hold on;
    plot(f,u,':m ');
    
    r=1;
    Fl=0;
    Fh=xlen/2;
    
    while r==1
    if filetype == 3   
        [X,Y,B] = ginput(1);
        XX(1)=X/Df;
        BB(1)=B;
        [X,Y,B]= ginput(1)
        XX(2)=X/Df;
        BB(2)=B;
        X=XX;
        B=BB(1);
        if floor(X(1))<((winlen+1)/2)
            X(1)=(winlen+1)/2;
        end
        if ceil(X(1))>(xlen-(winlen+1)/2)
            X(1)=xlen-( winlen+1)/2;
        end
        if floor(X(2))<((winlen+1)/2)
            X(2)=(winlen+1)/2;
        end
        if ceil(X(2))>(xlen-(winlen+1)/2)
            X(2)=xlen-( winlen+1)/2;
        end
        X
    else
        [X,Y,B] = ginput(1);
        X=X/Df;
        if floor(X)<((winlen+1)/2)
            X=(winlen+1)/2;
        end
        if ceil(X)>(xlen-(winlen+1)/2)
            X=xlen-( winlen+1)/2;
        end
    end
    if B~=1
        r=0;
    else
        if filetype==1
            %低通
            Fh=ceil(X);
            u=ones(1, xlen/2)*Amp;
            for i=1:winlen
                u(Fh-((winlen-1)/2+1)+i)=w(i);
            end
            for i=(Fh+(winlen-1)/2+1):xlen/2
                u(i)=0;
            end
        end
        
        if filetype==2            %高通
            Fl=ceil(X);
            u=ones(1, xlen/2)*Amp;
            for i=1:winlen
                u(Fl-((winlen-1)/2+1)+i)=v(i);
            end
            for i=1:(Fl-((winlen-1)/2+1))
                u(i)=0;
            end
        end
    
        if filetype==3            %带通
            temp=min(X);
            if temp>Fl
                 Fl=ceil(temp);
                 for i=1:winlen
                     u(Fl-((winlen-1)/2+1)+i)=v(i);
                 end
                 for i=1:(Fl-((winlen-1)/2+1))
                    u(i)=0;
                 end
            end
            temp=max(X);
            if Fh>temp
                Fh=ceil(temp);
                for i=1:winlen
                     u(Fh-((winlen-1)/2+1)+i)=w(i);
                end
                for i=(Fh+(winlen-1)/2+1):xlen/2
                    u(i)=0;
                end
            end
        end
    end
    clf;
    plot(f,z);
    title('信号的幅值谱图');
    xlabel('频率/Hz');
    ylabel('幅值/mv');
    hold on;
    plot(f,u,'m: ');
    end
    
    %% 计算滤波器频谱及脉冲响应函数
    disp('计算滤波器频谱及脉冲响应函数.....');
    u=u/Amp;
    s=fliplr(u);
    s=[u s];
    s=ifft(s);
    handle3 =figure;
    plot(real(s(1:100)));
    % N=input('输人滤波器单位脉冲响应序列的长度:');  % 如100
    N=100;  % 如100  输人滤波器单位脉冲响应序列的长度
    
    N=N/2;
    
    disp('进行傅里叶反变换,求单位脉冲响应函数');
    clf;
    s=fliplr(u);
    s=[u s];
    s=ifft(s);
    disps=[s(xlen-N:xlen) s(1:N)];
    plot(real(disps),'b')
    title('滤波器的单位脉冲响应函数');
    xlabel('采样点数');
    ylabel('幅值/mv');
    
    % disp('按任意键继续......');
    % pause
    
    %% 显示滤波结果
    disp('按任意键显示滤波结果......');
    handle4 =figure;
    subplot(2,1,1)
    plot(t, x, 'g');
    xlabel('时间/s');
    ylabel('原始信号');
    title('滤波结果比较');
    
    grid on;
    
    ft=fliplr(u);
    ft=[u, ft];
    ft=ft';
    ft= ft*(xlen/2);
    ft=absy.*ft;
    ft=ft.*exp(j*angley);
    fs= ifft(ft);
    subplot(2,1,2);
    plot(t,real(fs),'b');
    ylabel('滤波信号');
    xlabel('时间/s');
    grid on;
    % disp('按任意键保存滤波结果......');
    % pause
    
     %% 显示所有的窗口
    % figure(handle1);
    % figure(handle2);
    % figure(handle3);