阶次分析 程序 和 问题

sxf

信号：x = 0.7.*sin(2.*pi.*t.*t);
转速：n=t（转/秒）不是分钟
问题：
应该是 一阶 的
为什么我写的程序是：二阶的
谢谢！
调用方法：
clc;
clear;
L = 5000;                      % Length of signal
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                       % Sample time
t = (0:L-1)*T;
x = 0.7.*sin(2.*pi.*t.*t);
y=resample(x,0.005);   % 0.005 是重采样转数增量
L=length(y);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
Fs=1/0.005;
f = Fs/2*linspace(0,1,NFFT/2);
plot(f,2*abs(Y(1:NFFT/2)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
M文件resample:

function [y]=resample(signal,del)
L=length(signal);
j=0;
s=0.0;
while s<5.0
    s=sqrt(j*del*2);
    j=j+1;
end
N=j;
for i=1:N
    out(i)=sqrt(i*del*2);
end
for n=1:N
    for m=1:L-1
        if 0.001*m== out(n)
            ss(n)=signal(m);
            break;
        end
        if 0.001*m < out(n) && out(n)<0.001*(m+1)
           ss(n)=(signal(m+1)-signal(m))*(out(n)-0.001*m)/0.001+signal(m);
            break
        end
    end
end
y=ss;

sxf

今天想了想：是不是因为这是一个 调频信号  实际的得到信号的 频率是  2*t ？ 所以 实际就是 二阶 的 呢？

jkjxl

我是初学者，看不懂，楼主牛

loverpipi

一个正余弦信号的频率计算是对cos或者sin里面的多项式求导除以2*pi，如cos（2*pi*t*f）的频率就是f，而cos（2*pi*s（t））的频率就是对s（t）求导

紫砂清壶

%--------------------------------------------------------------------------
%-----数据导入(手动修改)--------------------------------------------
%--------------------------------------------------------------------------
N=2200000;
Fs=40000;
order_max=48;
t=(0:N-1)*1/Fs;
array_time_amp=Data1; %时域振动信号导入
pluse=Data2; %脉冲信号导入
%--------------------------------------------------------------------------
%-----计算脉冲发生时刻---------------------------------------------
%--------------------------------------------------------------------------
Num_pluse1=1;
threshold=200; %设定脉冲阈值
for temp2=1:length(pluse)-1;
if (pluse(temp2)<200&&pluse(temp2+1)>=200)
       Num_pluse1=[Num_pluse1,temp2+1];
end
end
if length(Num_pluse1)<2, return,end;
t_pluse=(Num_pluse1-1)/Fs;
%----------------------------------------------------------
%-------数字跟踪滤波-----------------------------------
%-----------------------------------------------------------
i=2;
while i<=length(t_pluse);
ft=2/(t_pluse(i)-t_pluse(i-1));
[b,a]=ellip(6,3,80,ft*order_max/(Fs/2));
array_time_amp(Num_pluse1(i-1):(Num_pluse1(i)-1))=filter(b,a,array_time_amp(Num_pluse1(i-1):(Num_pluse1(i)-1)));
i=i+1;
end
%----------------------------------------------------------
%---------重采样----------------------------------------
%----------------------------------------------------------
%------等角度时间计算--------------------------------
delt_thet=pi/24; %设定重采样角度
t_angle=[];
for temp3=3:length(t_pluse);
b=inv([1,t_pluse(temp3-2),t_pluse(temp3-2)^2;1,t_pluse(temp3-1),t_pluse(temp3-1)^2;1,t_pl
use(temp3),t_pluse(temp3)^2])*[0,2*pi,4*pi]';
if temp3==3;
k=0;
while k<1.5*2*pi/delt_thet;
tt=(sqrt(4*b(3)*(k*delt_thet-b(1))+b(2)^2)-b(2))/(2*(b(3)+eps));
t_angle=[t_angle,tt];
k=k+1;
end
else
k=pi/delt_thet;
while k>=pi/delt_thet && k<1.5*2*pi/delt_thet;
tt=(sqrt(4*b(3)*(k*delt_thet-b(1))+b(2)^2)-b(2))/(2*(b(3)+eps));
t_angle=[t_angle,tt];
k=k+1;
end
end
end
array_angle=[1:length(t_angle)].*delt_thet;
%------等角度时刻的幅值拟合（三种拟合方式）-----
%------线性插值拟合-----------------------------
% array_angle_amp=[];
% for temp5=1:length(t_angle);
% for temp4=1:length(t);
% if t_angle(temp5)>=t(temp4)&&t_angle(temp5)<t(temp4+1);
% %linear interpolation
% y=[array_time_amp(temp4),array_time_amp(temp4+1)];
% x=[t(temp4),t(temp4+1)];
% angle_amp=interp1(x,y,t_angle(temp5),'linear');
% array_angle_amp=[array_angle_amp,angle_amp];
%
% % -------- 三次多项式插值拟合------------------
% % if temp4==1;
% % y=[array_time_amp(temp4),array_time_amp(temp4+1)];
% % x=[t(temp4),t(temp4+1)];
% % angle_amp=interp1(x,y,t_angle(temp5),'linear');
% % array_angle_amp=[array_angle_amp,angle_amp];
% % else
% %
y=[array_time_amp(temp4-1),array_time_amp(temp4),array_time_amp(temp4+1),array_time_amp(temp4+2)];
% %由前两点与后两点得用三次多项式插值求得此角度对应的幅值
% % x=[t(temp4-1),t(temp4),t(temp4+1),t(temp4+2)];
% % angle_amp=interp1(x,y,t_angle(temp5),'cubic');
% % array_angle_amp=[array_angle_amp,angle_amp];
% % end
% break;
% end
% end
% end
%----------三次样条插值拟合-------------------------------
array_angle_amp=interp1(t,array_time_amp,t_angle,'cubic');
%--------------------------------------------------------
%------------------加窗--------------------------------
%-------------------------------------------------------
w=hanning(length(array_angle_amp))';
array_angle_amp=array_angle_amp.*w;
s%---------------------------------------------------------------
%-------------- FFT ------------------------------------------
%---------------------------------------------------------------
angle_dom_ffty=abs(fft(array_angle_amp))*2/length(array_angle_amp);
delt_order=2*pi/(length(angle_dom_ffty)*delt_thet);
angle_dom_fx=(0:length(angle_dom_ffty)-1)*delt_order;
FFTy=abs(fft(array_time_amp))*2/length(array_time_amp);
fx=(0:length(array_time_amp)-1)/t_sample;
%----------------------------------------------------------------
%---------------2 维图形显示----------------------------------
%----------------------------------------------------------------
figure(1);
subplot(2,1,1),plot(t,array_time_amp),title('time domianat'),xlabel('time'),ylabel('amplitude');
grid on;
subplot(2,1,2),plot(t,pluse),title('keyphasor'),xlabel('time'),ylabel('amplitude');
grid on;
figure(2);
subplot(2,1,1),plot(fx(1:length(fx)/2),FFTy(1:length(fx)/2)),title('Frequency domain');
xlabel('f'), ylabel('amplitude');
subplot(2,1,2),plot(array_angle,array_angle_amp), title('angle dominant'), xlabel('angle /rad'),
ylabel('amplitude');
grid on;
figure(3);
subplot(2,1,1),plot(angle_dom_fx(1:length(angle_dom_fx)/2),angle_dom_ffty(1:length(angle_dom_fx)/2));

impulse

紫砂清壶 发表于 2012-12-7 15:35

                               
登录/注册后可看大图

%--------------------------------------------------------------------------
%-----数据导入(手动修改 ...

这是一个离线分析程序，不能自动根据键相信号来计算阈值大小，而是采用人为设定的方法，还有键相信号的脉冲宽度也是确定的（为2）。

gogspe

2*t 是什么意思,请高手指导

plcplc

紫砂清壶 发表于 2012-12-7 15:35

                               
登录/注册后可看大图

%--------------------------------------------------------------------------
%-----数据导入(手动修改 ...

程序怎么不运行？我是未入门者，想学习。运行后出现：？？ Undefined function or variable 'Data1'.

gogspe

plcplc 发表于 2012-12-16 20:33

                               
登录/注册后可看大图

程序怎么不运行？我是未入门者，想学习。运行后出现：？？ Undefined function or variable 'Data1'.

Data1要自己输入

plcplc

gogspe 发表于 2012-12-17 21:27

                               
登录/注册后可看大图

Data1要自己输入

？whut？求qq。。。。

impulse

这个程序提供计算阶比跟踪算法的一个详细说明，应该说步骤很详细，但其处理效果欠佳。
下图是一个模拟信号的比较结果：
1、原始信号（频率从5Hz线性增加到30Hz，并且含有基频、0.5倍频、2倍频成分）

2、计算得到的转速曲线

3、该matlab算法得到的重采样波形与原始波形比较

4、采用NI的算法得到的重采样波形与原始波形比较

5、原始信号频谱与重采样信号阶比谱比较

6、附件为全部数据（float二进制数据），signal_x.dat为原始时域波形数据，signal_tacho.dat为键相脉冲数据，evenwave.dat为重采样后的波形数据，eventime.dat为重采样后的波形数据对应的时间序列。 data.rar (503.24 KB, 下载次数: 167)

2012-12-18 20:40 上传

点击文件名下载附件

firecat_2

impulse 牛啊

impulse

firecat_2 发表于 2013-1-6 10:44

                               
登录/注册后可看大图

impulse 牛啊

过奖，谢谢！我现在正在将这些算法和lms和bk的阶比跟踪算法就行比较。

impulse

计算阶比跟踪的难点其实在实时处理上，后处理在算法上不难，难的是实时，还要保证连续处理的波形是完全连续的，每一段波形前后不存在畸变。

westrongmc

impulse 发表于 2012-12-18 20:13

                               
登录/注册后可看大图

这个程序提供计算阶比跟踪算法的一个详细说明，应该说步骤很详细，但其处理效果欠佳。
下图是一个模拟信号 ...

有个疑问——为何0.5阶、2阶，有阶次模糊现象？
阶次分析是角度域内整周期采样，按理说，应该没有泄露，如上图中的第1阶所示。