求助 关于噪声信号的处理 声压 A声级 1/3倍频程
先把大致的思路讲一下1.获取仪器设备采集回来的噪声信号，这时的信号是电压信号2.将上述电压信号除以仪器的灵敏度获取声压信号3.将上述的声压信号进行fft快速傅里叶变换，频率的上限为采样频率的一半(这里的fft变换采用intel fortran中mkl里的库函数)4.对获取的fft结果用abs()函数求绝对值，再把绝对值乘以2除以数据的个数n获取幅值5.根据国家标准给定的标称中心频率计算1/3倍频程的上下限频率，确定各个频带6.分别对每个频带进行累加求和，求和的方式是幅值的平方求和再除以27.再对每个频带求和结果进行对数运算，10*log(频带求和结果/p0^2)，其中p0=2*10^(-5)8.再对对数运算后的各个频带的结果进行A计权加权，加权值参照国家标准中给出的权值9.导出上诉运算结果
由于程序是用fortran编写的，就不把程序贴出来了，下面把对 1000Hz 94dB 灵敏度为51.8mV/Pa的噪声源采集并处理后的结果贴出来10.0000000000000 -102.029494712842
12.5000000000000 -108.425002105719
16.0000000000000 -92.2024358413517
20.0000000000000 -94.8378861356631
25.0000000000000 -85.3236931670592
31.5000000000000 -72.0851717728487
40.0000000000000 -69.1864488051986
50.0000000000000 -51.5691890371075
63.0000000000000 -76.5416584798834
80.0000000000000 -74.5817738068465
100.000000000000 -67.4763276711592
125.000000000000 -48.3164556745454
160.000000000000 -39.1097124346171
200.000000000000 -39.7814848232638
250.000000000000 -51.1406939098685
315.000000000000 -29.4861408303704
400.000000000000 -49.5664368499122
500.000000000000 -48.4382820296058
630.000000000000 -37.1939040211018
800.000000000000 -25.4736445671135
1000.00000000000 76.5281679598751
1250.00000000000 -20.9315652339948
1600.00000000000 -30.2028169383888
2000.00000000000 -10.2249031076519
2500.00000000000 -32.2811288452300
3150.00000000000 -19.5323774166911
4000.00000000000 -29.5006162356157
5000.00000000000 -29.8800717995616 从上述计算结果可以看出结果存在明显的偏差，1000Hz时的结果仅为76.5281679598751
希望各位版友能结合前面提供的思路帮助看下问题出在哪里有哪里没看明白的可以再详细解释这里先谢谢各位版友 没有人回应么？ lz很用心学习，先鼓励一个。看了下你的处理过程，思路上没什么问题，差异这么大很可能是你程序写的有错误吧，另外，建议不要用纯音验证你的算法过程，1000Hz纯音没法验证你的倍频程计算和A计权计算过程，用宽带谱噪声吧。
上程序，我帮你看看。 回复 4 # hyl2323 的帖子
先谢版主了 program fspg
use MKL_DFTI
integer i,j,n1,n2
complex,allocatable ::sp(:),tt(:),xx(:,:)
double precision,allocatable ::amp(:),ff(:),fst(:),freal(:),aw(:),spg(:),spl(:)
double precision temp,p0,fs,finala,sensitivity
type(DFTI_DESCRIPTOR), POINTER :: My_Desc1_Handle
integer :: Status
open(unit=1,file='data.dat',status='old',access='sequential',form='formatted')
open(unit=2,file='awei.dat',status='old',access='sequential',form='formatted')
open(unit=3,file='aweighting.dat',status='old',access='sequential',form='formatted')
write(*,*)'请输入采集数据总数'
read(*,*)n1
write(*,*)'请输入灵敏度(mV/Pa)'
read(*,*)sensitivity
allocate(xx(n1,2),sp(n1),tt(n1),amp(n1),ff(n1),fst(n2),freal(n2*2),aw(n2),spg(n2),spl(n1))
do 100,i=1,n1,1
read(1,*)xx(i,1),xx(i,2)
tt(i)=xx(i,2)/sensitivity
sp(i)=0.0
100 continue
fs=n1/(xx(n1,1)-xx(1,1))
write(*,*)'信号采样频率为',fs,'Hz'
do 110,i=1,34,1
read(2,*)fst(i),aw(i)
110 continue
n2=0
do 120,i=1,34,1
if(fs/2.0 .gt. fst(i))then
n2=n2+1
endif
120 continue
temp=0.0
p0=2.0*10**(-5.0)
do 130,i=1,n1,1
temp=temp+tt(i)
130 continue
do 140,i=1,n1,1
tt(i)=tt(i)-temp/n1
140 continue
Status=DftiCreateDescriptor(My_Desc1_Handle,DFTI_SINGLE,DFTI_COMPLEX,1,n1)
Status=DftiSetValue(My_Desc1_Handle,DFTI_PLACEMENT,DFTI_NOT_INPLACE)
Status=DftiCommitDescriptor(My_Desc1_Handle)
Status=DftiComputeForward(My_Desc1_Handle,tt,sp)
Status=DftiFreeDescriptor(My_Desc1_Handle)
do 150,i=1,n1/2,1
amp(i)=abs(sp(i))*2.0/n1
ff(i)=(i-1.0)*fs/n1
150 continue
temp=0.0
do 160,j=1,n2,1
freal(2*j-1)=(1.0/2.0**(1.0/6.0))*fst(j)
freal(2*j)=2.0**(1.0/6.0)*fst(j)
do 170,i=1,n1/2,1
if(freal(2*j-1) .lt. ff(i) .and. freal(2*j) .gt. ff(i))then
temp=temp+amp(i)**2/2.0
endif
170 continue
write(*,*)temp/p0**2
spg(j)=10*log(temp/p0**2)+aw(j)
temp=0.0
160 continue
do 180,i=1,n2,1
write(*,*)fst(i),spg(i)
write(3,*)fst(i),spg(i)
180 continue
temp=0.0
do 190,i=1,n2,1
temp=temp+10.0**(0.1*(spg(i)))
190 continue
finala=10.0*log(temp)
write(*,*)'合成A声级'
write(*,*)finala
deallocate(xx,sp,tt,amp,ff,fst,freal,aw,spg,spl)
close(1)
close(2)
close(3)
pause
end program fspg
open(unit=1,file='data.dat',status='old',access='sequential',form='formatted')
打开设备采集的数据文件
open(unit=2,file='awei.dat',status='old',access='sequential',form='formatted')
打开a声级计权频率和权值文件
open(unit=3,file='aweighting.dat',status='old',access='sequential',form='formatted')
打开最后保存数据的文件
do 100,i=1,n1,1
read(1,*)xx(i,1),xx(i,2)
tt(i)=xx(i,2)/sensitivity
sp(i)=0.0
100 continue
读入采集的数据文件并除以灵敏度
do 110,i=1,34,1
read(2,*)fst(i),aw(i)
110 continue
读入a声级计权频率和权值
do 130,i=1,n1,1
temp=temp+tt(i)
130 continue
do 140,i=1,n1,1
tt(i)=tt(i)-temp/n1
140 continue
消去直流分量
Status=DftiCreateDescriptor(My_Desc1_Handle,DFTI_SINGLE,DFTI_COMPLEX,1,n1)
Status=DftiSetValue(My_Desc1_Handle,DFTI_PLACEMENT,DFTI_NOT_INPLACE)
Status=DftiCommitDescriptor(My_Desc1_Handle)
Status=DftiComputeForward(My_Desc1_Handle,tt,sp)
Status=DftiFreeDescriptor(My_Desc1_Handle)
快速傅里叶变换
do 150,i=1,n1/2,1
amp(i)=abs(sp(i))*2.0/n1
ff(i)=(i-1.0)*fs/n1
150 continue
修正幅值并计算对应频率
do 160,j=1,n2,1
freal(2*j-1)=(1.0/2.0**(1.0/6.0))*fst(j)
freal(2*j)=2.0**(1.0/6.0)*fst(j)
do 170,i=1,n1/2,1
if(freal(2*j-1) .lt. ff(i) .and. freal(2*j) .gt. ff(i))then
temp=temp+amp(i)**2/2.0
endif
170 continue
write(*,*)temp/p0**2
spg(j)=10*log(temp/p0**2)+aw(j)
temp=0.0
160 continue
分频带累加并进行权值相加
写出了部分注释，版友可以参考下
这个帖子希望大家多多讨论
再次感谢版主
纠正已经发现了的错误
log()应改为log10()
这个是我疏忽了……………… 本帖最后由 hyl2323 于 2011-3-7 19:02 编辑
log修正为log10后，那结果不是更小了啊？
灵敏度有没有化mV为V？
不能跟你调试，看不到中间结果，发现不了错误，还得靠你自己一步步验证中间结果吧。
上两张计算完以后用结果画的图
版友们看看有没有什么问题
hyl2323 发表于 2011-3-7 18:48 static/image/common/back.gif
log修正为log10后，那结果不是更小了啊？
灵敏度有没有化mV为V？
不能跟你调试，看不到中间结果，发现不了 ...
确实……采集的时候电压是V，后来老师告诉我是mV，纠结了好长时间终于弄出来了
我已经上了后来计算的结果的图，版主再帮忙看看，谢啦 恭喜lz成功了。{:{23}:} 恭喜楼主成功了。
本来想好好学习，但fortan语言没有怎么用过，看的有点糊涂 楼主，后辈无知也遇到了同样的问题，看了你的编程思路，对第四步中“再把绝对值乘以2除以数据的个数n获取幅值”和第6步中“再除以2”不明白，去年的这个问题你早已解决了吧，还望不吝赐教。 回复 1 # coolzlj 的帖子
楼主，后辈无知也遇到了同样的问题，看了你的编程思路，对第四步中“再把绝对值乘以2除以数据的个数n获取幅值”和第6步中“再除以2”不明白，去年的这个问题你早已解决了吧，还望不吝赐教。
页:
[1]
2