求助 关于噪声信号的处理 声压 A声级 1/3倍频程

coolzlj

先把大致的思路讲一下

1.获取仪器设备采集回来的噪声信号，这时的信号是电压信号

2.将上述电压信号除以仪器的灵敏度获取声压信号

3.将上述的声压信号进行fft快速傅里叶变换，频率的上限为采样频率的一半

(这里的fft变换采用intel fortran中mkl里的库函数)

4.对获取的fft结果用abs()函数求绝对值，再把绝对值乘以2除以数据的个数n获取幅值

5.根据国家标准给定的标称中心频率计算1/3倍频程的上下限频率，确定各个频带

6.分别对每个频带进行累加求和，求和的方式是幅值的平方求和再除以2

7.再对每个频带求和结果进行对数运算，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     

希望各位版友能结合前面提供的思路帮助看下问题出在哪里  

有哪里没看明白的可以再详细解释  

这里先谢谢各位版友

coolzlj

没有人回应么？

hyl2323

lz很用心学习，先鼓励一个。看了下你的处理过程，思路上没什么问题，差异这么大很可能是你程序写的有错误吧，另外，建议不要用纯音验证你的算法过程，1000Hz纯音没法验证你的倍频程计算和A计权计算过程，用宽带谱噪声吧。

hyl2323

上程序，我帮你看看。

coolzlj

回复 4 # hyl2323 的帖子

先谢版主了

coolzlj

    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

coolzlj

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
分频带累加并进行权值相加


写出了部分注释，版友可以参考下
这个帖子希望大家多多讨论
再次感谢版主

coolzlj

纠正已经发现了的错误
log()应改为log10()
这个是我疏忽了………………

hyl2323

log修正为log10后，那结果不是更小了啊？
灵敏度有没有化mV为V？
不能跟你调试，看不到中间结果，发现不了错误，还得靠你自己一步步验证中间结果吧。

coolzlj

上两张计算完以后用结果画的图
版友们看看有没有什么问题

coolzlj

hyl2323 发表于 2011-3-7 18:48

                               
登录/注册后可看大图

log修正为log10后，那结果不是更小了啊？
灵敏度有没有化mV为V？
不能跟你调试，看不到中间结果，发现不了 ...

确实……采集的时候电压是V，后来老师告诉我是mV，纠结了好长时间终于弄出来了
我已经上了后来计算的结果的图，版主再帮忙看看，谢啦

hyl2323

恭喜lz成功了。

3QMM

恭喜楼主成功了。
本来想好好学习，但fortan语言没有怎么用过，看的有点糊涂

ahjcs24_2006

楼主，后辈无知也遇到了同样的问题，看了你的编程思路，对第四步中“再把绝对值乘以2除以数据的个数n获取幅值”和第6步中“再除以2”不明白，去年的这个问题你早已解决了吧，还望不吝赐教。

ahjcs24_2006

回复 1 # coolzlj 的帖子

楼主，后辈无知也遇到了同样的问题，看了你的编程思路，对第四步中“再把绝对值乘以2除以数据的个数n获取幅值”和第6步中“再除以2”不明白，去年的这个问题你早已解决了吧，还望不吝赐教。