声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1515|回复: 7

[综合讨论] 请帮我看看这个频谱分析问题做的对吗?

[复制链接]
发表于 2007-8-4 22:38 | 显示全部楼层 |阅读模式

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

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

x
我对一个地震波首先进行了去噪分析,又进行了速度-位移fft变换和频谱分析。请各位帮我看看我做的对吗?另外,作阈值去噪分析和滤波的区别是什么?如果做0.2~50Hz的带通滤波,该怎么做啊?

load JGS-UD(40-45);
[swa,swd]=swt(x,2,'db1');
[thr,sorh]=ddencmp('den','wv',x);
dswd=wthresh(swd,sorh,thr);
clean=iswt(swa,dswd,'db1');
it=1; sf=100; n=length(clean);
t=0:1/sf:(n-1)/sf; %建立时间向量
nfft=2^nextpow2(n); %大于并最接近n的2的幂次方为FFT长度
y=fft(clean,nfft); %FFT变换
m=20*log10(abs(y)*2/n);
ff=(0:nfft/2)*sf/nfft;
df=sf/nfft; %计算频率间隔(Hz/s)
dw=2*pi*df; %计算圆频率间隔(rad/s)
w1=0:dw:2*pi*0.5*sf; %建立正的离散圆频率向量
w2=-2*pi*(0.5*sf-df):dw:-dw; %建立负的离散圆频率向量
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
y=ifft(a,nfft); %IFFT变换
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n)); Y=fft(y,nfft); M=2*abs(Y)/n; Pyy=M.^2;
figure(1);
subplot(211),plot(t,clean),title('De-noised signal'); xlabel('时间'),ylabel('速度');grid on;
subplot(212),plot(t,y),title('位移图'); xlabel('时间'),ylabel('位移');grid on;
figure(2); loglog(ff,M(1:257)),title('幅值频谱图'); xlabel('频率');grid on;
figure(3); loglog(ff,Pyy(1:257)),title('能量频谱图'); xlabel('频率'); grid on;

[ 本帖最后由 ChaChing 于 2010-3-6 14:18 编辑 ]

JGS-UD(40-45).mat

2.2 KB, 下载次数: 22

回复
分享到:

使用道具 举报

 楼主| 发表于 2007-8-5 14:09 | 显示全部楼层

补充图像

补充图像
1.jpg
2.jpg
3.jpg
发表于 2007-8-5 18:44 | 显示全部楼层
看了一下,程序是对的。
但对要做0.2~50Hz的带通滤波有些疑问。从上谱图可看出最高频率为50Hz,带通上限设在50Hz没有意义;又在谱分析中ff(2)是接近0.2Hz,把带通下限设在0.2Hz同样没有多大的意义,0.2~50Hz的带通滤波变成了全通。同时在谱图上看不到低于0.2Hz有明显干扰,为什么楼主要设置这样的带通呢?
 楼主| 发表于 2007-8-5 19:20 | 显示全部楼层
那之前我发的那个帖子中为了消除低频振荡而采用的0.2~50Hz的带通滤波,单是把原程序加入
fmin=0.2; fmax=50Hz;
%计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
a=zeros(1,nfft);
%消除指定正频带外的频率成分
a(ni:na)=y(ni:na);
%消除指定负频带外的频率成分
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
我把这些加到原程序里了,但好像位移和速度的图像是一样的,没得到您所做出 的结果?原帖地址:
http://forum.vibunion.com/forum/thread-48861-1-1.html
发表于 2007-8-5 20:35 | 显示全部楼层
这次的情形和上次并不一样,这次数据只有512个点,所以得到我上帖提到的分辨率,ff(2)是接近0.2Hz。而在http://forum.vibunion.com/forum/thread-48861-1-1.html帖子上是处理JGS-UD.txt数据,数据长10000,在FFT处理时nfft=16384,这样分辨率就完全不同,ff(2)=0.006Hz,0.2Hz对应的是笫33条谱线,在这样情况下进行滤波是有意义的,而在ff(2)是接近0.2Hz下要把带通下限设在0.2Hz是没有意义的。所以对于具体问题要作具体的分析。

评分

1

查看全部评分

 楼主| 发表于 2007-8-5 21:39 | 显示全部楼层
在上帖中您为了消除低频振荡,采用了0.2~50Hz的带通滤波,我做了一下,可是为什么位移的时程曲线和速度时程曲线一模一样呢?频谱图倒是有变化。
x=load('JGS-UD.txt');
x=x';
[swa,swd]=swt(x,3,'db1');
[thr,sorh]=ddencmp('den','wv',x);
dswd=wthresh(swd,sorh,thr);
clean=iswt(swa,dswd,'db1');
fmin=0.2;
fmax=50;
it=1;
sf=100;
n=length(clean);
%建立时间向量
t=0:1/sf:(n-1)/sf;
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%FFT变换
y=fft(clean,nfft);
ff=(0:nfft/2)*sf/nfft;
%计算频率间隔(Hz/s)
df=sf/nfft;
%计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
%计算圆频率间隔(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
a=zeros(1,nfft);
%消除指定正频带外的频率成分
a(ni:na)=y(ni:na);
%消除指定负频带外的频率成分
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
%IFFT变换
y=ifft(a,nfft);
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n));
Y=fft(y,nfft);
M=2*abs(Y)/n;
Pyy=M.^2;
figure(1);
subplot(211),plot(t,clean),title('De-noised signal');
xlabel('时间'),ylabel('速度');grid on;
subplot(212),plot(t,y),title('位移图');
xlabel('时间'),ylabel('位移');grid on;
figure(2);
loglog(ff,M(1:8193)),title('幅值频谱图');
xlabel('频率');grid on;
figure(3);
plot(ff,Pyy(1:8193)),title('能量频谱图');
xlabel('频率');grid on;
又是哪里出现问题了呢?找了很久也没找出原因,又得麻烦您了。

JGS-UD.txt

38.04 KB, 下载次数: 14

发表于 2007-8-6 07:51 | 显示全部楼层

a=a1-a2*i;
改为
y=a1-a2*i;
便可以了。
 楼主| 发表于 2007-8-6 08:06 | 显示全部楼层
谢谢
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-29 00:04 , Processed in 0.074162 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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