%?這是時域積分的程序，原始信號為速度信號y1
%?x_velocity??速度信號
%?D?數據長度
%?time_series?分析時間序列
%?total_time??采樣時間
%mean_velocity??速度均值
%?spec_area?速度頻譜
%freq_area?分析頻率序列

%有效值?Xrms
close?all
clear?all
clc;
%=========================================================================%
?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%??信號產生???%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%?
%=========================================================================%
N=1024;
fs=250;
fc=40;
fb=80;
t=（0:1/fs:（N-1）/fs）;
y0=10*sin（2*pi*fc*t+200*pi/180）+20*cos（2*pi*fb*t+150*pi/180）+6*sin（2*pi*0.5*t+150*pi/180）+0.15;??%模擬信號作為加速度信號
%加白噪聲，a=0;b=0.5;??%均值為a，方差為b^2
%y2=a+b*randn（1length（y1））?高斯白噪聲?title（‘N（0，0.25）的高斯白噪聲‘）;
n=length（y0）;
y1=wgn（1n4）;?????????????????????????????????????????????????????????????%產生高斯白噪聲
y=y0+y1;???????????????????????????????????????????????????????????????????%加入噪聲之后的信號
x_accr=y;
D=length（x_accr）;??????????????????????????????????????????????????????????%數據長度
time_series=（1/fs:1/fs:D/fs）;??????????????????????????????????????????????%分析時間序列也可time_series=（0:D-1）/fs;????????????????????????????????????????????????%分析時間序列的另一種表示
time_series=time_series‘;?????????????????????????????????????????
total_time=D/fs;???????????????????????????????????????????????????????????%采樣時間，作為橫坐標的范圍
N=2^nextpow2（D）;???????????????????????????????????????????????????????????%確定分析頻率范圍的點數:取最接近D的較大2次冪

spec_y0=fft（y0N）;
spec_y=fft（yN）;
freq_series=fs*（0:N/2-1）/N;

figure（4）
subplot（311）;plot（ty1‘r‘）;
axis（[0?1.1?-5??5]）;
ylabel（‘y（um）‘）;xlabel（‘t（s）‘）;
title（‘強度為4dwb的高斯白噪聲‘）;
grid?on;
figure（1）
subplot（421）;plot（ty0）;
ylabel（‘y0（um）‘）;xlabel（‘t（s）‘）;
title（‘理想信號‘）;
grid?on;
subplot（422）;plot（freq_seriesabs（spec_y0（1:N/2））*2/N）;
xlabel（‘頻率（Hz）‘）;
ylabel（‘幅值（mm/g）‘）;
title（‘理想加速度幅頻譜‘）;
grid?on;
subplot（423）;plot（ty1）;
ylabel（‘y（um）‘）;xlabel（‘t（s）‘）;
title（‘原始信號（含噪）‘）;
grid?on;
subplot（424）;plot（freq_seriesabs（spec_y（1:N/2））*2/N）;
xlabel（‘頻率（Hz）‘）;
ylabel（‘幅值（mm/g）‘）;
title（‘原始加速度幅頻譜‘）;
grid?on;

%=========================================================================%
?%%%%%%%%%%%%%%%%%%%%%%%%%??加速度信號濾波處理???%%%%%%%%%%%%%%%%%%%%%%%%%%%
%=========================================================================%?????????????????????????????????????????????????????????
mean_accr=mean（x_accr）;????????????????????????????????????????????????????%去直流
disp（‘實際的直流量‘）;?????????????????????????
disp（mean_accr）;
x_accr=x_accr-mean_accr;????
spec_accr=fft（x_accrN）;???????????????????????????????????
freq_series=fs*（0:N/2-1）/N;????????????????????????????????????????????????%建立分析頻率序列
%?輸出加速度時域波形和頻譜圖
subplot（425）;
plot（time_seriesx_accr‘r‘）;
ax=max（abs（x_accr））;
ra=ax+1;
axis（[0?total_time?-ra?ra]）;
xlabel（‘采樣時間（s）‘）;
ylabel（‘幅值（mm/g）‘）;
title（‘振動加速度，已去直流‘）;????????????????????
grid?on;
subplot（426）;
plot（freq_seriesabs（spec_accr（1:N/2））*2/N‘r‘）;?
xlabel（‘頻率（Hz）‘）;
ylabel（‘幅值（mm/g）‘）;
title（‘去直流的加速度幅頻譜‘）;
grid?on;
%設計高通濾波器