本帖最后由 牛小贱 于 2015-3-8 16:14 编辑
我改了圆频率向量之后，用标准的正弦信号加上一个常数项，得出的位移有上升的趋势，下面是我的程序：
- function PY1(xx)
- %频域积分法求速度及位移
- fs=1000; %采样频率
- % ff=20;
- % N=1000;
- % tt=(0:N-1)/fs;
- % xx=3*sin(2*pi*ff*tt)+0.02;
- fmin=0.2; %最小截止频率
- fmax=500; %最大截止频率
- c=9.88e+003; %单位变换系数，位移单位为mm
- it=2; %积分次数
- n=length(xx);
- t=0:1./fs:(n-1)./fs; %建立时间向量
- nfft=2^nextpow2(n); %FFT变换长度，大于并最接近n的2的幂次方为FFT的长度
- y=fft(xx,nfft); %FFT变换
- subplot(3,1,2); %绘制积分前的时程曲线图形
- plot(t,xx,'k-');
- xlabel('时间（s）');
- ylabel('加速度幅值（g）');
- title('图2 加速度时程曲线图');
- grid on;
- %=========================================================================%
- %%%%%%%%%%%%%%%%%%%%% 绘制频谱图，计算机床振动频率 %%%%%%%%%%%%%%%%%%%%
- %=========================================================================%
- P=y.*conj(y)/nfft; %求取功率密度，其中conj(y)是求y的共轭
- f=fs*(0:(nfft/2-1))/nfft; %设定频率变化范围
- subplot(3,1,1);
- plot(f,P(1:nfft/2),'b-'); %绘制频谱图
- xlabel('频率（Hz）');
- ylabel('功率密度');
- title('图1 功率密度频谱图');
- grid on;
- %=========================================================================%
- %%%%%%%%%%%%%%%%%%%%%% 进行积分的频域变换 %%%%%%%%%%%%%%%%%%%%%
- %=========================================================================%
- df=fs/nfft; %计算频率间隔（Hz/s）
- ni=round(fmin/df+1); %计算指定频带对应频率数组的下标
- na=round(fmax/df+1);
- dw=2*pi*df; %计算圆频率间隔（rad/s）
- w1=0:dw:2*pi*0.5*fs; %建立正的离散圆频率向量
- w2=-2*pi*(0.5*fs-df):dw:-dw; %建立负的离散圆频率向量
- % w1=0:dw:2*pi*(0.5*fs-df);
- % w2=2*pi*(0.5*fs-df):-dw:0;
- w=[w1,w2]; %将正负圆频率向量组合成一个向量
- w=w.^it; %以积分次数为指数，建立圆频率变量向量
- a=zeros(1,nfft); %进行积分的频域变换
- a(2:(nfft-1))=y(2:(nfft-1))./w(2:(nfft-1));
- % if it==2
- % y=-a; %进行二次积分的相位变换，二次积分因为是-w^2，所以有y=-a
- % else %进行一次积分的相位变换，一次积分时考虑到jw，把实部和虚部对换
- % a1=imag(a);
- % a2=real(a);
- % y=a1-a2*i;
- % end
- y=-a; %进行二次积分的相位变换
- a=zeros(1,nfft);
- a(ni:na)=y(ni:na); %消除指定正频带外的频率成分
- a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1); %消除指定负频带外的频率成分
- y=ifft(a,nfft); %IFFT变换
- y=real(y(1:n))*c; %取逆变换的实部n个元素并乘以单位变换系数为积分结果
- %%
- % subplot(3,1,2); %绘制积分前的时程曲线图形
- % plot(t,xx,'k-');
- % xlabel('时间（s）');
- % ylabel('加速度幅值（g）');
- % title('图2 加速度时程曲线图');
- % grid on;
- subplot(3,1,3); %绘制积分后的时程曲线图形
- plot(t,y,'r-');
- xlabel('时间（s）');
- ylabel('位移幅值（mm）');
- title('图3 位移时程曲线图');
- grid on;
- end
|
[复制链接]
发表于 2014-6-18 15:37
显身卡
