FSWT与TT.rar_FSWT与TT_matlab FSWT下载_故障诊断_轴承故障诊断_频率切片小波故障诊断

%FSWT-TT对实际震动信号的处理： close all clear all [t1,h]=textread('H:\故障信号\近电机端轴承外圈故障\第一次测试\采样频率12000\转速1450rpm\恒转速\测试(2017-03-03,15-48-57)\DataRecor3.txt','%f%f','headerlines',38); h=h(1:8192); N=length(h); y=h; fs=12000;%采样频率在此 t=0:1/fs:(N-1)/fs;% 将以下程序copy到matlab编辑器中运行，或直接在工作区运行即可 SNR=5; y=awgn(y,SNR, 'measured', [], 'dB'); load('D:\信号处理桌面集合\FSWT-TT变换论文相关\TT变换程序中的数据\real1.mat') y=h; figure() plot(t,y); xlim([0,0.68]); % axis([0,inf,-4,5]) xlabel('Time(s)') ylabel('A/(m.s^-^2)') title('轴承故障仿真信号时域波形图') xlabel('Time(s)') ylabel('A/(m.s^-^2)') s=y; N=length(s); Title1={'原始信号的频谱'}; MyFFT(s,N,fs,1,Title1);%figure H1=abs(hilbert(h)); H=H1-mean(H1); b=abs(fft(H)); f=fs*(0:N-1)/N; b=b*2/N; figure(); plot(f,b); xlim([0,1000]); xlabel('f/Hz'); ylabel('A/(m.s^-^2)'); title('原始信号包络结果');6 f1=0; f2=1200; %%% [f1,f2] is your observed frequency range, you can change any scope k1=fix(f1*N/fs-0.5); k2=fix(f2*N/fs-0.5); df=1; %% is a observed frequency step length if(k2>N/2+1) k2=N/2+1; end%K2不能超过信号长度的一半 fp= fix(k1:df:k2); %% fp if the observed frequency in discrete form kapa=sqrt(2)/2/0.025; %%kapa is the time-frequency resolution factor Tn=N; %% smaple numer in time domain, you can chage it [A] = GetFSWT(s,fs,fp,kapa,Tn); %SET_Y_exact ss1 = GetInvFSWT(N,A,fp);%recover the signal from the band of fp; f1=1201; f2=2400; %%% [f1,f2] is your observed frequency range, you can change any scope k1=fix(f1*N/fs-0.5); k2=fix(f2*N/fs-0.5); df=1; %% is a observed frequency step length if(k2>N/2+1) k2=N/2+1; end%K2不能超过信号长度的一半 fp= fix(k1:df:k2); %% fp if the observed frequency in discrete form kapa=sqrt(2)/2/0.025; %%kapa is the time-frequency resolution factor Tn=N; %% smaple numer in time domain, you can chage it [A] = GetFSWT(s,fs,fp,kapa,Tn); %SET_Y_exact ss2 = GetInvFSWT(N,A,fp);%recover the signal from the band of fp; f1=2401; f2=3600; %%% [f1,f2] is your observed frequency range, you can change any scope k1=fix(f1*N/fs-0.5); k2=fix(f2*N/fs-0.5); df=1; %% is a observed frequency step length if(k2>N/2+1) k2=N/2+1; end%K2不能超过信号长度的一半 fp= fix(k1:df:k2); %% fp if the observed frequency in discrete form kapa=sqrt(2)/2/0.025; %%kapa is the time-frequency resolution factor Tn=N; %% smaple numer in time domain, you can chage it [A] = GetFSWT(s,fs,fp,kapa,Tn); %SET_Y_exact ss3 = GetInvFSWT(N,A,fp);%recover the signal from the band of fp; f1=3601; f2=4800; %%% [f1,f2] is your observed frequency range, you can change any scope k1=fix(f1*N/fs-0.5); k2=fix(f2*N/fs-0.5); df=1; %% is a observed frequency step length if(k2>N/2+1) k2=N/2+1; end%K2不能超过信号长度的一半 fp= fix(k1:df:k2); %% fp if the observed frequency in discrete form kapa=sqrt(2)/2/0.025; %%kapa is the time-frequency resolution factor Tn=N; %% smaple numer in time domain, you can chage it [A] = GetFSWT(s,fs,fp,kapa,Tn); %SET_Y_exact ss4 = GetInvFSWT(N,A,fp);%recover the signal from the band of fp; f1=4800; f2=fs/2; %%% [f1,f2] is your observed frequency range, you can change any scope k1=fix(f1*N/fs-0.5); k2=fix(f2*N/fs-0.5); df=1; %% is a observed frequency step length if(k2>N/2+1) k2=N/2+1; end%K2不能超过信号长度的一半 fp= fix(k1:df:k2); %% fp if the observed frequency in discrete form kapa=sqrt(2)/2/0.025; %%kapa is the time-frequency resolution factor Tn=N; %% smaple numer in time domain, you can chage it [A] = GetFSWT(s,fs,fp,kapa,Tn); %SET_Y_exact ss5 = GetInvFSWT(N,A,fp);%recover the signal from the band of fp; figure() subplot(511) plot((0:N-1)/fs,ss1); % ylabel('A/(m.s^-^2)'); % xlabel('t/s'); xlim([0,0.68]); title={''}; hold on subplot(512) plot((0:N-1)/fs,ss2); % ylabel('A/(m.s^-^2)'); % xlabel('t/s'); xlim([0,0.68]); title={''}; hold on subplot(513) plot((0:N-1)/fs,ss3); ylabel('A/(m.s^-^2)'); % xlabel('t/s'); xlim([0,0.68]); title={''}; hold on subplot(514) plot((0:N-1)/fs,ss4); % ylabel('A/(m.s^-^2)'); % xlabel('t/s'); xlim([0,0.68]); title={''}; hold on subplot(515) plot((0:N-1)/fs,ss5); % ylabel('A/(m.s^-^2)'); xlabel('t/s'); xlim([0,0.68]); title={''}; hold on % % % % % B=sqrt(A.*conj(A)); % % % % % mx= max(max(B)); % % % % % B=fix(B*128/mx); % % % % % t= (0:Tn-1)*N/fs/Tn; % figure() % plot((0:N-1)/fs,ss1,'r'); KI1=kurt(ss1,'kurt2'); KI2=kurt(ss2,'kurt2'); KI3=kurt(ss3,'kurt2'); KI4=kurt(ss4,'kurt2'); KI5=kurt(ss5,'kurt2'); % % % % % figure() % % % % % image(fp*fs/N,t,B); %%%%%二维图 % % % % % figure(); % % % % % [x,y]=meshgrid(fp*fs/N,t); % % % % % [na,nb]=size(x); % % % % % z=reshape(B,na,nb); % % % % % mesh(x,y,z); % % % % % C=zeros(Tn,Tn); % % % % % for i=1:Tn % % % % % C(i,:)=ifft(A(i,:),N);%%%%%%%%%%% % % % % % end % % % % % D=real(C );%%%%%% % % % % % figure() % % % % % mesh(D); % % % % % figure() % % % % % image(Tn,t,D); % % % % % D=fliplr(D); % % % % % TTreshape=diag(D); % % % % % % % % % % % plot(t,TTreshape); % % % % % figure() % % % % % plot(t,TTreshape); % % % % % xlim([0,0.68]); % % % % % ylim([-2000,2000]); % % % % % xlabel('Time(s)') % % % % % ylabel('Amplitude') % % % % % title('TT变化对角线提取的时域结果'); % % % % % % % % % % Title1={'FSWT-tt后对角线信号的频谱'}; % % % % % MyFFT(TTreshape,N,fs,1,Title1); % % % % % % % % % % H1=abs(hilbert(TTreshape)); % % % % % H=H1-mean(H1); % % % % % b=abs(fft(H)); % % % % % f=fs*(0:N-1)/N; % % % % % b=b*2/N; % % % % % figure();%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % plot(f,b);%降噪后重构波形的Hilbert包络谱 % % % % % xlim([0,1000]); % % % % % xlabel('f/Hz'); % % % % % ylabel('A/(m.s^-^2)'); % % % % % title('对角线直接包络结果'); % % % % % % % % % % % ssssnr=snr(ss1); % % % % % H1=abs(hilbert(ss1)); % % % % % H=H1-mean(H1); % % % % % b=abs(fft(H)); % % % % % f=fs*(0:N-1)/N; % % % % % b=b*2/N; % % % % % figure();%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % plot(f,b);%降噪后重构波形的Hilbert包络谱 % % % % % xlim([0,1000]); % % % % % xlabel('f/Hz'); % % % % % ylabel('A/(m.s^-^2)'); % % % % % title('选频包络结果'); % % % % % % % % % % % KI=kurt(ss1,'kurt2'); % %%降噪部分：： % [U,S,V]=svd(D); %求R的奇异值,同时保存求奇异值所需的两个酉矩阵 % A = diag(S); %保存所有奇异值,并保存奇异值的个数 % figure() %输出奇异值序列 % plot(A(1:200)); % xlabel('奇异值序号Sn'); % ylabel('奇异值V'); % N1=Tn;N2=Tn; % NewS=zeros(N1,N2); % for i=1:1:11 % % NewS(i,:)=S(i,:); % end % %对保留的奇异值，重构信号 % NewMatrix=U*NewS*V'; % % TTreshape_jiangzao=diag(NewMatrix); % figure() % plot(t,real(TTreshape_jiangzao)) % title('TT变换后 矩阵降噪信号的频谱 奇异值=32的时域'); % Title1={'TT变换后 矩阵降噪信号的频谱 奇异值=32的频谱'}; % MyFFT(TTreshape_jiangzao,N,fs,1,Title1); % % % H1=abs(hilbert(TTreshape)); % H=H1-mean(H1); % b=abs(fft(H)); % f=fs*(0:N-1)/N; % b=b*2/N; % figure();%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % plot(f,b);%降噪后重构波形的Hilbert包络谱 % xlim([0,400]); % xlabel('f/Hz'); % ylabel('A/(m.s^-^2)'); % title('对角线直接包络结果'); % % % % H1=abs(hilbert(TTreshape)); % H=H1-mean(H1); % b=abs(fft(H)); % f=fs*(0:N-1)/N; % b=b*2/N; % figure();%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % plot(f,b);%降噪后重构波形的Hilbert包络谱 % xlim([0,400]); % xlabel('f/Hz'); % ylabel('A/(m.s^-^2)'); % title('对角线直接包络结果'); % % % % % % %对对角线直接降噪 % %% SVD处理 对对角线直接降噪。 % xtn=TTreshape;%TTreshape; % k=length(xtn); %重构矩阵维数 % Data=zeros(1,k); % H=myhankel01(xtn); % [U,S,V]=svd(H); %求R的奇异值,同时保存求奇异值所需的两个酉矩阵 % A = diag(S); %保存所有奇异值,并