matlab-实现变转速信号的阶次分析和变转速信号角域重采样

%% Order Analysis_1(the method is differ to the paper) % Order Signal Fs =8000; t = (0:160000)/Fs; Ticks = 10; % Number of Impulses per revolution of light barrier rpstheo = (t.^6.*exp(-t/1))+10; % Theoretical rotational speed in rev. per sec. angle4ticks = cumsum(2*pi*rpstheo)/Fs; % Total angle = Integrate(ang. veloc.)|0...t encodersig = (1+sin(angle4ticks*Ticks))>=1; % Calculation of encoder signal % Create a measurement signal with some fix, variing and noisy parts: backnoise = 0.1.*conv(randn(1,length(encodersig)-11),[-1 0 2 -3 5 4 0 0 0 -2 0 0]); backsignal = 0.1*sin(t).*sin(2*pi*2*pi*13*t)+... 0.2*cos(t.^1.5).*cos(2*pi*2*pi*31*t)+... +0.15*sin(2*pi*2*pi*59*t); ordersignal = 1*sin(t*2.9).*sin(angle4ticks*1)+... 0.5*sin(t*2.9).*sin(angle4ticks*3)+... 0.3*sin(t*2.9).*sin(angle4ticks*5)+... 0.1*sin(t*2.9).*sin(angle4ticks*7)+... 0.4*abs(cos(angle4ticks/5)).*sin(angle4ticks*11.2); signal = backnoise+backsignal+ordersignal; signal = detrend(signal);% 消除线性趋势. F_signal = fft(signal); f=(1:80000)/8000; figure plot(abs(F_signal(1:80000))/8000); % Signal Reconstruction rps = rpstheo; phi = cumsum(2*pi*rps/Fs); Nphin = length(phi)*1; phin = linspace(0.1,phi(end),Nphin)'; % Interpolate signaln = interp1(phi,signal,phin,'linear'); F_signaln=fft(signaln(2:160001)); figure plot(f,abs(F_signaln(1:80000))/8000); figure subplot(211);plot(abs(F_signal(1:80000))/8000); subplot(212);plot(f,abs(F_signaln(1:80000))/8000); %% Order Analysis_2(the method is differ to the paper) % Order Signal % clc; % clear; Fs =8000; t = (0:160000)/Fs; rpstheo = 4.*t; % Theoretical rotational speed in rev. per sec. %注意当rpstheo =t+a,a较t的值越大,角域重采样效果越差. angle4ticks = cumsum(2*pi*rpstheo)/Fs; % Total angle = Integrate(ang. veloc.)|0...t % Create a measurement signal with some fix, variing and noisy parts: backnoise =5.* randn(size(t)); ordersignal = sin(2*pi*10.*(4.*t.^2)); signal =backnoise+ordersignal; signal = detrend(signal);% 消除线性趋势. F_signal = fft(signal); f =(1:80000)/400; figure plot(f,abs(F_signal(1:80000))/800); % Signal Reconstruction rps = rpstheo; phi = cumsum(2*pi*rps/Fs); Nphin = length(phi); phin = linspace(0.1,phi(end),Nphin)'; % Interpolate signaln = interp1(phi,signal,phin,'spline'); F_signaln =fft(signaln); figure plot(f,abs(F_signaln(1:80000))/800); % PSD psd=(abs(F_signaln(1:80000))/800).^2/20; figure plot(f,psd); subplot(511);plot(t,signal); subplot(512);plot(f,abs(F_signal(1:80000))/800); subplot(513);plot(t,signaln); subplot(514);plot(f,abs(F_signaln(1:80000))/800); subplot(515);plot(f,psd); %axis([10 100 -0.5 0.5]); %set(gcf,'units','centimeters','position',[6 6 14 5]); %% 阶比仿真 重采样角度间隔不理想,且效果仍然不理想,有待改进(提升角域分辨率). % clear all % close all % clc Fs=256;Ts=1/Fs;N=2048;t=0*Ts:Ts:(N-1)*Ts; %采样时间为8秒 x=sin(2*pi*(t.^2))+sin(4*pi*(t.^2))+sin(6*pi*(t.^2))+sin(10.2*pi*(t.^2))+0.8*randn(size(t)); y1=(abs(fft(x)))*2/N;f1=(0:N/2-1)*Fs/N; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fc=128;%理想低通滤波器截止频率 wc=2*fc*pi;%256*pi %转速信号，1024转 %n=2*60*t.^2;%转速 i=1; a=zeros(32);%等时间的时刻值序列 for nn=1:32 % if nn==1 % a(1)=0; %else a(nn)=sqrt(nn*2); %end % i=i+1; end b=[0;2*pi;4*pi]; kk=1; y=zeros(1,12);%等角度时刻序列 for i=1:2:length(a)-3 A=[1,a(i),a(i)^2;1,a(i+1),a(i+1)^2;1,a(i+2),a(i+2)^2]; % b1=[2*pi;4*pi;6*pi]; b2=inv(A)*b; for k=100:299 %采样角度大于等于pi小于等于3*pi y(kk)=(sqrt(4*b2(3)*(k*pi/100-b2(1))+b2(2)^2)-b2(2))/b2(3)/2; kk=kk+1; end end yy=zeros(1,length(y));%重构信号名称 for ii=1:length(y) SR=0; for m=1:N SR=SR+x(m)*Ts*sin(wc*(y(ii)-m*Ts))/(y(ii)-m*Ts); %x(m)为函数所有时刻对应值的遍历,m*Ts为所有时刻遍历 end yy(ii)=SR/pi;% 上述循环完成低通滤波器的功能 end %阶比横坐标计算 N2=length(yy); fs=100; ff=(0:N2/2-1)*fs/N2; y3=abs(fft(yy))*2/N; figure subplot(411);plot(f1,y1(1:N/2)); subplot(412);plot(ff,y3(1:N2/2)); subplot(413);plot(x); subplot(414);plot(yy) %% 阶比跟踪算法 较前面两种方法使用更方便,但结果尚不理想. %xtn为输出：等角度采样的信号序列 %输入：x为等时间间隔采样信号序列，t为时间，fs采样频率，Dmax为最大阶次，pf为转速曲线序列,order:拟合转速曲线的阶次，wu：舍弃的点数 %function [Tn,xtn] = getCOT(x,t,fs,Dmax,pf,order,wu) % clear all % close all % clc Fs=256;Ts=1/Fs;N=2048;t=0*Ts:Ts:(N-1)*Ts; %采样时间为8秒 x=sin(2*pi*(t.^2))+sin(4*pi*(t.^2))+sin(6*pi*(t.^2))+sin(10.2*pi*(t.^2))+0.8*randn(size(t)); y=(abs(fft(x)))*2/N;f1=(0:N/2-1)*Fs/N; fc=128;%理想低通滤波器截止频率 wc=2*fc*pi;%256*pi figure plot(t,x); figure plot(abs(fft(x))); Dmax=1000; pf=t; order=1; wu=10; t = t - min(t); dw = pi/Dmax; %重采样角度间隔 dt = 1/Fs ;%采样时间间隔 a = polyfit(t,pf,order); %3阶拟合：ft = a(1)*t.^3 + a(2)*t.^2 + a(3)*t+a(4); ft = polyval(a,t); %得到拟合后的频率曲线序列 Na = length(a); for j = 1:Na a(j) = a(j)/(Na-j+1);%对转速函数进行积分,得到的函数系数. end lenXtn = fix(2*pi*sum(ft*dt)/dw); %计算重采样后的数据长度,sum(ft*dt)表示对转过角度的累计量,即总旋转角度量. lenXtn = lenXtn-wu; Tn = zeros(1,lenXtn); % 计算键相时标 ,即重采样角域等间隔坐标. for ii = 1 :lenXtn % 求解方程 temp = ii/(2*Dmax); r = roots([a -temp]); for kk = 1: length(r) if isreal(r(kk)) if r(kk) > 0 && r(kk) < 10000 Tn(ii) = r(kk); end end end end xtn=zeros(1,lenXtn); for ii=1:length(Tn) SR=0; for m=1:2048 SR=SR+x(m)*Ts*sin(256*pi*(Tn(ii)-m*Ts))/(Tn(ii)-m*Ts); %x(m)为函数所有时刻对应值的遍历,m*Ts为所有时刻遍历 end xtn(ii)=SR; % 上述循环完成低通滤波器的功能 end figure plot(Tn,xtn); figure plot(abs(fft(xtn))); axis([1 10000 0 10^5.1]); figure subplot(411);plot(t,x); subplot(412);plot(abs(fft(x))); subplot(413);plot(Tn,xtn); subplot(414);plot(abs(fft(xtn)));