请各位帮我分析一下作出来的频谱图好像不对啊？（上传了数据）

鸭鸭

我在作一个地震波的频谱分析，该数据为时域的速度数据。我将其转成了时域的位移，然后作频谱分析，请大家帮我分析一下为什么频谱图不对？数据 BYS.dat (14.75 KB, 下载次数: 175)

2007-7-20 16:12 上传

点击文件名下载附件

或者 BYS.txt (14.75 KB, 下载次数: 152)

2007-7-20 16:12 上传

点击文件名下载附件

。
程序如下：

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');

复制代码

songzy41

从楼主的程序中没有看到哪一语句是在求频谱图的。从量纲上来看，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和它的频谱图。

鸭鸭

首先谢谢sonzy41回复。还有些疑问：
那有没有必要将速度转成位移，再来求位移的频谱呢？如果要做位移频谱如何做？为什么将x和y的坐标轴分别定成0-50和-60-20呢？有什么依据吗

鸭鸭

> 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；
这样对吗？为什么位移频谱和速度频谱几乎是一样的 啊？

songzy41

原帖由 鸭鸭 于 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（单位的分贝）是为了图中最大幅图地显示图形。

鸭鸭

看了帖子http://forum.vibunion.com/forum/thread-47704-1-1.html中给出的从加速度到速度，到位移的程序，从加速度到位移的变换，仅仅是将程序中的it改成2吗？我试了一下，好像改不改得出的都是加速度和速度的图像啊，得不到位移的。请您指点！还有sx、sy1、sy2的标注是怎么弄上去的，也没有另外输入啊？是先将时域的加速度作积分变换成频域的速度，然后再逆变换得到速度的时域曲线，基本思想是这样么？我初学有很多地方不懂，希望SONZY41不吝赐教！:loveliness:

鸭鸭

还有您给的程序最后输出积分后的数据时，数据仍保存在同样的文件名下，就把原数据给替换了，想请教一下将该程序怎样修改，才能既输出积分后数据又能保留原始数据资料呢？

鸭鸭

%建立正的离散圆频率向量
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帮我编写一下位移的频谱分析程序吧，我学习一下

songzy41

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

songzy41

原帖由 鸭鸭 于 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就可以了。

鸭鸭

%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;
我做了一下，但感觉好像不对啊，请老师批改作业！

songzy41

主要程序中y=a1-a2*i和y=ifft(a,nfft)不一致造成的错误。改成a=a1-a2*i以后就好了。在笫3幅图上横坐标最好用频率，而不是圆频率。

鸭鸭

%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;
谢谢谢谢，谢谢您的耐心。这样对了么？频谱的线怎么那么密啊

junzifei

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

[ 本帖最后由 junzifei 于 2007-7-28 23:00 编辑 ]

songzy41

我在http://forum.vibunion.com/forum/viewthread.php?tid=47704&highlight=帖子中提到了楼主最后计算的结果位移中有一低频振荡，其原因也已在该帖中说明了。因为楼主是新手，所以在程序中让他把滤波部分简化掉了。如果还是按我在该帖中讲到的增加一个低频滤波，在本帖子中设为0.2~50Hz的带通滤波，这样便能得到如下图的数据，位移的低频振荡完全消除了。