滤波程序.m

matlab编写滤波的源代码

clear;
clc;
fni=input('频域积分-输入数据文件名：','s');
load(fni);
fs=20000;               %采样频率
n=24000;                 %数据点数
n1=24000;
n2=48000;
A=1000;                   %放大倍数为1mv/m/s2
dt =1/fs;
t1=Data1_time_AI_0(1:n);%信号1的时间轴
x1=Data1_AI_0(1:n);     %信号1
Wp = 1000/10000; %低通滤波
Ws =1050/10000;
Rp=3;
Rs = 60;
[N, Wn] = buttord(Wp, Ws, Rp, Rs);
n3=6;
[b,a] = butter(n3,Wn);
y1=filter(b,a,x1);
rms1=sqrt(mean((y1*A).*(y1*A)));%信号rms值
L1=20*log10(rms1/0.000001);%信号振动加速度级
peak1=sort(y1*A,'descend');
figure(1);
axis fill;
plot(t1,y1*A);
xlabel('Time (s)');
ylabel('Amplitude (m/s^2)');
axis fill;
title('3.5Mpa 测点CJ9 y方向 加速度波形图(1000HZ以下)');
grid on;
xlim([0 1.2]);

t2=Data1_time_AI_0(1:n);%信号1的时间轴
x2=Data1_AI_0(1:n);     %信号1  ％带通滤波1000－3000HZ
Wp = [1000 3000]/10000;
Ws = [950 3050]/10000;
Rp=3;
Rs = 60;
[N, Wn] = buttord(Wp, Ws, Rp, Rs);
n3=6;
[b,a] = butter(n3,Wn);
y2=filter(b,a,x2);
rms2=sqrt(mean((y2*A).*(y2*A)));%信号rms值
L2=20*log10(rms2/0.000001);%信号振动加速度级
peak2=sort(y2*A,'descend');
figure(2);
plot(t2,y2*A);
xlabel('Time (s)');
ylabel('Amplitude (m/s^2)');
axis fill;
title('3.5Mpa 测点CJ9 y方向 加速度波形图(1000~3000Hz)');
grid on;
xlim([0 1.2]);

t3=Data1_time_AI_0(1:n);%信号1的时间轴
x3=Data1_AI_0(1:n);     %信号1  ％带通滤波3000－6000HZ
Wp = [3000 6000]/10000;
Ws = [2950 6050]/10000;
Rp=3;
Rs = 60;
[N, Wn] = buttord(Wp, Ws, Rp, Rs);
n3=6;
[b,a] = butter(n3,Wn);
y3=filter(b,a,x3);
rms3=sqrt(mean((y3*A).*(y3*A)));%信号rms值
L3=20*log10(rms3/0.000001);%信号振动加速度级
peak3=sort(y3*A,'descend');
figure(3);
plot(t3,y3*A);
xlabel('Time (s)');
ylabel('Amplitude (m/s^2)');
axis fill;
title('3.5Mpa 测点CJ9 y方向 加速度波形图(3000~6000Hz)');
grid on;
xlim([0 1.2]);

t4=Data1_time_AI_0(1:n);%信号1的时间轴
x4=Data1_AI_0(1:n);     %信号1  ％高通滤波6000HZ以上
Wp = 6000/10000;
Ws =5950/10000;
Rp=3;
Rs = 60;
[N, Wn] = buttord(Wp, Ws, Rp, Rs);
n3=6;
[b,a] = butter(n3,Wn,'high');
y4=filter(b,a,x4);
rms4=sqrt(mean((y4*A).*(y4*A)));%信号rms值
L4=20*log10(rms4/0.000001);%信号振动加速度级
peak4=sort(y4*A,'descend');
figure(4);
plot(t4,y4*A);
xlabel('Time (s)');
ylabel('Amplitude (m/s^2)');
axis fill;
title('3.5Mpa 测点CJ9 y方向 加速度波形图(6000Hz以上)');
grid on;
xlim([0 1.2]);

b1=rms1/(rms1+rms2+rms3+rms4);
b2=rms2/(rms1+rms2+rms3+rms4);
b3=rms3/(rms1+rms2+rms3+rms4);
b4=rms4/(rms1+rms2+rms3+rms4);

%带阻滤波
t1=Data1_time_AI_1(1:10000);%信号1的时间轴
x1=Data1_AI_1(n1:n);     %信号1
Wp = [950 3050]/5000;
Ws = [1000 3000]/5000;
Rp=3;
Rs = 13;
[N, Wn] = BUTTORD(Wp, Ws, Rp, Rs);
n3=6;
[b,a] = butter(n3,Wn,'stop');%带阻滤波器
y1=filter(b,a,x1);