声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 10679|回复: 26

[FFT] 请各位帮我分析一下作出来的频谱图好像不对啊?(上传了数据)

[复制链接]
发表于 2007-7-20 16:12 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
本帖最后由 牛小贱 于 2014-6-29 13:20 编辑

我在作一个地震波的频谱分析,该数据为时域的速度数据。我将其转成了时域的位移,然后作频谱分析,请大家帮我分析一下为什么频谱图不对?数据 BYS.dat (14.75 KB, 下载次数: 175) 或者 BYS.txt (14.75 KB, 下载次数: 152)
程序如下:
  1. load BYS.dat;
  2. v=BYS';
  3. t=0:0.01:39.99;
  4. y=fft(v);
  5. N=length(y);
  6. f=100*(1:N)/N;
  7. s=y./(j*(2*pi*f));
  8. Y=ifft(s);
  9. S=real(Y);
  10. m=abs(Y);
  11. Pyy=abs(Y(1:N)).^2;
  12. subplot(4,1,1),plot(v,'r');
  13. subplot(4,1,2),plot(S,'g');
  14. subplot(4,1,3),loglog(f,m,'m');
  15. subplot(4,1,4),plot(f,m,'k');
复制代码


回复
分享到:

使用道具 举报

发表于 2007-7-20 16:58 | 显示全部楼层
从楼主的程序中没有看到哪一语句是在求频谱图的。从量纲上来看,v是在时间域的,y和s是在频率域的,Y又是在时间域的,S和m都是从Y而得,都是时间域的信号。y是频率域,可由它来求出信号v的频谱。
load BYS.dat;
v=BYS';
t=0:0.01:39.99;
y=fft(v);
N=length(y);
f=100*(1:N)/N;
subplot(2,1,1),plot(t,v,'r'); grid;
subplot 212; plot(f,20*log10(abs(y)*2/N));
axis([0 50 -60 20]);  grid;
给出了信号v和它的频谱图。
 楼主| 发表于 2007-7-20 17:37 | 显示全部楼层
首先谢谢sonzy41回复。还有些疑问:
那有没有必要将速度转成位移,再来求位移的频谱呢?如果要做位移频谱如何做?为什么将x和y的坐标轴分别定成0-50和-60-20呢?有什么依据吗
 楼主| 发表于 2007-7-20 17:44 | 显示全部楼层
> load BYS.dat;
v=BYS';
t=0:0.01:39.99;
y=fft(v);
N=length(y);
f=100*(1:N)/N;
s=y./(j*(2*pi*f));
subplot(2,1,1),plot(t,v,'r'); grid;
subplot 212; plot(f,20*log10(abs(s)*2/N));
axis([0 60 -100 0]);  grid;
这样对吗?为什么位移频谱和速度频谱几乎是一样的 啊?
发表于 2007-7-20 20:22 | 显示全部楼层
本帖最后由 wdhd 于 2016-6-3 10:04 编辑
原帖由 鸭鸭 于 2007-7-20 17:37 发表
首先谢谢sonzy41回复。还有些疑问:
那有没有必要将速度转成位移,再来求位移的频谱呢?如果要做位移频谱如何做?为什么将x和y的坐标轴分别定成0-50和-60-20呢?有什么依据吗

从速度到位移需要积分,在帖子http://www.chinavib.com/forum/thread-47704-1-1.html中给出了从加速度到速度,到位移的程序,楼主可参考一下该程序,修改变成速度到位移的程序。
楼主在程序中给了t,间隔为0.01,表示采样频率为100,因此在求频谱时最高只能到50,这便是频谱中为什么x轴取0-50的原因;y轴取-60-20(单位的分贝)是为了图中最大幅图地显示图形。
 楼主| 发表于 2007-7-23 09:19 | 显示全部楼层
看了帖子http://forum.vibunion.com/forum/thread-47704-1-1.html中给出的从加速度到速度,到位移的程序,从加速度到位移的变换,仅仅是将程序中的it改成2吗?我试了一下,好像改不改得出的都是加速度和速度的图像啊,得不到位移的。请您指点!还有sx、sy1、sy2的标注是怎么弄上去的,也没有另外输入啊?是先将时域的加速度作积分变换成频域的速度,然后再逆变换得到速度的时域曲线,基本思想是这样么?我初学有很多地方不懂,希望SONZY41不吝赐教!:loveliness:
 楼主| 发表于 2007-7-23 10:06 | 显示全部楼层
还有您给的程序最后输出积分后的数据时,数据仍保存在同样的文件名下,就把原数据给替换了,想请教一下将该程序怎样修改,才能既输出积分后数据又能保留原始数据资料呢?
 楼主| 发表于 2007-7-23 10:32 | 显示全部楼层
%建立正的离散圆频率向量
w1=0:dw:2*pi*0.5*sf;
%建立负的离散圆频率向量
w2=-2*pi*(0.5*sf-df):dw:-dw;
%将正负圆频率向量组合成一个向量
w=[w1,w2];
在速度积分成位移时为什么还要用到离散圆频率呢?普通的f=100*(1:N)/N;或nquist频率不行吗?离散圆频率又是什么?w1=0:dw:2*pi*0.5*sf中的0.5是怎么得来的?如何知道最大和最小截止频率?那个x和y的标注我明白是怎么弄的了,麻烦帮我解答一下这几个问题吧!谢谢
如果您觉得麻烦,或者您用matlab帮我编写一下位移的频谱分析程序吧,我学习一下
发表于 2007-7-23 11:23 | 显示全部楼层
原帖由 鸭鸭 于 2007-7-23 09:19 发表
看了帖子http://www.chinavib.com/forum/thread-47704-1-1.html中给出的从加速度到速度,到位移的程序,从加速度到位移的变换,仅仅是将程序中的it改成2吗?我试了一下,好像改不改得出的都是加速度和速度的图像啊,得不到位移的。请您指点!还有sx、sy1、sy2的标注是怎么弄上去的,也没有另外输入啊?是先将时域的加速度作积分变换成频域的速度,然后再逆变换得到速度的时域曲线,基本思想是这样么?我初学有很多地方不懂,希望SONZY41不吝赐教!

在帖子http://www.chinavib.com/forum/thread-47704-1-1.htmlunknowno是从王济的书中引用来的,建议楼主看一下该书。从加速度到速度的变换时设置it=1,从加速度到位移的变换时设置it=2。而且计算的结果速度和位移的图象是不一样的,在该帖的笫6层和笫8层看到幅值就差2个数量级。关于sx、sy1、sy2的标注楼主开始时可以不设置,在最后作图时也不去用这些参数。因为原程序是引用书中的,有些不符合我们实际的使用,所以上帖中我建议楼主在该程序基础上改一下。

评分

2

查看全部评分

发表于 2007-7-23 11:41 | 显示全部楼层
原帖由 鸭鸭 于 2007-7-23 10:06 发表
还有您给的程序最后输出积分后的数据时,数据仍保存在同样的文件名下,就把原数据给替换了,想请教一下将该程序怎样修改,才能既输出积分后数据又能保留原始数据资料呢?
原帖由 鸭鸭
在速度积分成位移时为什么还要用到离散圆频率呢?普通的f=100*(1:N)/N;或nquist频率不行吗?离散圆频率又是什么?w1=0:dw:2*pi*0.5*sf中的0.5是怎么得来的?如何知道最大和最小截止频率?那个x和y的标注我明白是怎么弄的了,麻烦帮我解答一下这几个问题吧!谢谢
如果您觉得麻烦,或者您用matlab帮我编写一下位移的频谱分析程序吧,我学习一下

在原程序中,读数据和写数据是在两个不同的文件中,楼主在程序中可看到以下两个语句,分别表示打开读数据和写数据的文件,是两个文件,所以没有把原始数据冲了。
fid=fopen(fni,'r');
fid=fopen(fno,'w');

又在原程序中用到圆频率w(=ω)=2*pi*f,所谓离散圆频率是因为f是离散值。因为有一个因子1/jω。如果用f也一样的,f要进行归一化,并乘2*pi,最后还是得到圆频率。
最大和最小截止频率这部分也可删除,原程序中是用来滤波的。
楼主的关键是由速度求出位移。知道位移后求谱值,只要用FFT就可以了。



点评

赞成: 5.0
赞成: 5
  发表于 2014-6-29 13:21

评分

1

查看全部评分

 楼主| 发表于 2007-7-23 16:26 | 显示全部楼层
%fni=input('频域积分-输入数据文件名:','s');
fni='JGS-UD.txt';
fid=fopen(fni,'r');
sf=fscanf(fid,'%f',1);%采样频率
c=fscanf(fid,'%f',1);%单位变换系数
it=fscanf(fid,'%f',1);%积分次数
sx=fscanf(fid,'%s',1);%横向坐标轴的标注
sy1=fscanf(fid,'%s',1);%纵向坐标轴输入单位的标注
sy2=fscanf(fid,'%s',1);%纵向坐标轴输出单位的标注
fno=fscanf(fid,'%s',1);%输出数据文件名
x=fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
it=1;
n=length(x);
%建立时间向量
t=0:1/sf:(n-1)/sf;
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%FFT变换
y=fft(x,nfft);
ff=(0:nfft/2)*sf/nfft;
%计算频率间隔(Hz/s)
df=sf/nfft;
%计算圆频率间隔(rad/s)
dw=2*pi*df;
%建立正的离散圆频率向量
%w1=0:dw:2*pi*0.5*sf;
w1=0:dw:2*pi*0.5*sf;
%建立负的离散圆频率向量
w2=-2*pi*(0.5*sf-df):dw:-dw;
%w2=2*pi*(0.5*sf-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;
else
%进行一次积分的相位变换
a1=imag(a);
a2=real(a);
y=a1-a2*i;
end
%IFFT变换
y=ifft(a,nfft);
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n))*c;
Y=fft(y,nfft);
M=20*log10(abs(Y)*2/n);

subplot(3,1,1);
plot(t,x);
xlabel(sx);
ylabel(sy1);
grid on;
subplot(3,1,2);
plot(t,y);
xlabel(sx);
ylabel(sy2);
grid on;
subplot(3,1,3); plot(w,M);
grid;
我做了一下,但感觉好像不对啊,请老师批改作业!
JGS-UD.jpg

JGS-UD.txt

38.09 KB, 下载次数: 86

发表于 2007-7-23 18:51 | 显示全部楼层
主要程序中y=a1-a2*i和y=ifft(a,nfft)不一致造成的错误。改成a=a1-a2*i以后就好了。在笫3幅图上横坐标最好用频率,而不是圆频率。

点评

赞成: 5.0
赞成: 5
  发表于 2014-6-29 13:21
回复 支持 1 反对 0

使用道具 举报

 楼主| 发表于 2007-7-23 20:36 | 显示全部楼层
%fni=input('频域积分-输入数据文件名:','s');
fni='JGS-UD.txt';
fid=fopen(fni,'r');
sf=fscanf(fid,'%f',1);%采样频率
c=fscanf(fid,'%f',1);%单位变换系数
it=fscanf(fid,'%f',1);%积分次数
sx1=fscanf(fid,'%s',1);%横向坐标轴的标注
sx2=fscanf(fid,'%s',1);%横向坐标轴的标注
sy1=fscanf(fid,'%s',1);%纵向坐标轴输入单位的标注
sy2=fscanf(fid,'%s',1);%纵向坐标轴输出单位的标注
sy3=fscanf(fid,'%s',1);%纵向坐标轴输出单位的标注
fno=fscanf(fid,'%s',1);%输出数据文件名
x=fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
it=1;
n=length(x);
%建立时间向量
t=0:1/sf:(n-1)/sf;
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%FFT变换
y=fft(x,nfft);
ff=(0:nfft/2)*sf/nfft;
%计算频率间隔(Hz/s)
df=sf/nfft;
%计算圆频率间隔(rad/s)
dw=2*pi*df;
%建立正的离散圆频率向量
%w1=0:dw:2*pi*0.5*sf;
w1=0:dw:2*pi*0.5*sf;
%建立负的离散圆频率向量
w2=-2*pi*(0.5*sf-df):dw:-dw;
%w2=2*pi*(0.5*sf-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;
else
%进行一次积分的相位变换
a1=imag(a);
a2=real(a);
a=a1-a2*i;
end
%IFFT变换
y=ifft(a,nfft);
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n))*c;
Y=fft(y,nfft);
M=20*log10(abs(Y)*2/n);

subplot(3,1,1);
plot(t,x);
xlabel(sx1);
ylabel(sy1);
grid on;
subplot(3,1,2);
plot(t,y);
xlabel(sx1);
ylabel(sy2);
grid on;
subplot(3,1,3); plot(ff,M(1:8193));
xlabel(sx2);
ylabel(sy3);
grid;
谢谢谢谢,谢谢您的耐心。这样对了么?频谱的线怎么那么密啊
JGS-UD2.jpg
发表于 2007-7-28 22:56 | 显示全部楼层

回复 #13 鸭鸭 的帖子

楼主原始数据去除趋势项了吗?
原始数据就是速度吗?
起始时很小速度的值代表什么?是测量噪声还是速度值。
积分后有上升趋势呀!

[ 本帖最后由 junzifei 于 2007-7-28 23:00 编辑 ]
发表于 2007-7-29 09:58 | 显示全部楼层
我在http://forum.vibunion.com/forum/viewthread.php?tid=47704&highlight=帖子中提到了楼主最后计算的结果位移中有一低频振荡,其原因也已在该帖中说明了。因为楼主是新手,所以在程序中让他把滤波部分简化掉了。如果还是按我在该帖中讲到的增加一个低频滤波,在本帖子中设为0.2~50Hz的带通滤波,这样便能得到如下图的数据,位移的低频振荡完全消除了。
ayy3c.jpg

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-23 21:23 , Processed in 0.080668 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表