[已解决] 瞬时声压时域数据怎么用matlab进行1/3倍频程声压级分析

  [复制链接]
王猫猫040 发表于 2015-6-15 11:04:00
LZ你好!我有几个问题想请教你
if fu>fs/2
         m=j-1;
         break;

这里的m是哪儿来的呢?
nc=30,但是我用你的程序执行了一下之后得到的length(yc)=13,这是怎么回事呢?正确的length(yc)应该等于30才对吧
麻烦LZ解答,谢谢!

王猫猫040 发表于 2015-6-16 14:56:01
simplebinbin 发表于 2015-3-23 09:12
方法1与声级计计算结果差不多

LZ你好!我有几个问题想请教你
if fu>fs/2
          m=j-1;
          break;

这里的m是哪儿来的呢?
nc=30,但是我用你的程序执行了一下之后得到的length(yc)=13,这是怎么回事呢?正确的length(yc)应该等于30才对吧
麻烦LZ解答,谢谢!

simplebinbin 发表于 2015-7-2 16:27:55
length(yc)=13是因为你的采样频率比较低,只能计算到中心频率为315Hz的频段
m用来计数,等于length(yc),绘图时可用,bar(Lp1(1;m));

simplebinbin 发表于 2015-7-2 16:28:48
王猫猫040 发表于 2015-6-16 14:56
LZ你好!我有几个问题想请教你
if fu>fs/2
          m=j-1;

length(yc)=13是因为你的采样频率比较低,只能计算到中心频率为315Hz的频段
m用来计数,等于length(yc),绘图时可用,bar(Lp1(1;m));

王猫猫040 发表于 2015-7-12 09:25:57
simplebinbin 发表于 2015-7-2 16:27
length(yc)=13是因为你的采样频率比较低,只能计算到中心频率为315Hz的频段
m用来计数,等于length(yc),绘 ...

谢谢你详细的回复~

我还有几个问题需要请教你:

1)由瞬时声压时程计算总声压级时,加汉宁窗这部分
w=hanning(n);
yy=1.633*y.*w;是什么意思呢?(或者说1.633这么奇怪的数字是哪里来的)

2)计算三分之一倍频对应的声压级时,
   b=zeros(1,nfft);
    b(nl:nu)=a(nl:nu);
    b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
    c=ifft(b,nfft);
    %计算对应每个中心频段的声压有效值
    yc(j)=sqrt(var(real(c(1:n))));
end
这几句程序字面上的意思我能明白,但是其物理意义是什么呢?

没学过信号处理,但是做课题需要用到这些程序,自己研究了好久,也没找到介绍如何由声压时程计算总声压级和A计权声压级的资料。。。表示很头疼。。麻烦卤煮解答~~

simplebinbin 发表于 2015-7-14 16:20:50
王猫猫040 发表于 2015-7-12 09:25
谢谢你详细的回复~

我还有几个问题需要请教你:

1、1.633是能量修正系数,加窗以后保证能量不变,网上搜一下,有相关资料
2、将1/3倍频程中心频率对应频带的值提取出来,通过IFFT变换,得到这个频带的时域声压值(相当于一个带通滤波器),在这里直接用带通滤波器也可以,滤波上下限为:nu/nl

王猫猫040 发表于 2015-7-17 16:39:14
simplebinbin 发表于 2015-7-14 16:20
1、1.633是能量修正系数,加窗以后保证能量不变,网上搜一下,有相关资料
2、将1/3倍频程中心频率对应频 ...

多谢!!!

buaasuozi 发表于 2015-10-20 22:17:07
simplebinbin 发表于 2014-5-21 10:39
参考别的文献,发现有新的方法计算噪声1/3倍频程计权声级的方法,参考文献:张登攀《噪声1/3倍频程计权声级 ...

非常感谢楼主分享资源
关于这个程序,是不是没有绘制完全?缺少A计权的

另外,我想问一下,f f1 fc cf都是怎么获取的呢?
有什么依据吗?

simplebinbin 发表于 2015-10-27 14:05:48
buaasuozi 发表于 2015-10-20 22:17
非常感谢楼主分享资源
关于这个程序,是不是没有绘制完全?缺少A计权的

1、未计权的数据出来了,加上A计权的值就得到计权了的
2、f,f1,fc为1/3之一倍频程中心频率,cf为A计权值参考标准HJ453-2008,IEC61672-1等

hebutkappa 发表于 2015-11-9 16:14:26
本帖最后由 hebutkappa 于 2015-11-9 06:30 编辑
songzy41 发表于 2014-5-14 03:40
LZ在声压计中对于p0肯定不是一个设置的值:p0=2e-5;%参考声压,而是通过标准信号发生器对了传声器发出一个 ...


请问宋老师:
我有一个频域-声压曲线。横坐标HZ纵坐标db
可否直接用下面程序转成1/3倍频啊~做出来的数值不对~~能不能指点下啊
f=[1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00];
fc=[f,10*f,100*f,1000*f,10000*f];
%中心频率与下限频率的比值
oc6=2^(1/6);
nc=length(fc);
n=length(x);
nfft=2^nextpow2(n);
a=fft(x,nfft);
for j=1:nc
    fl=fc(j)/oc6;
    fu=fc(j)*oc6;
    nl=round(fl*nfft/sf+1);
    nu=round(fu*nfft/sf+1);
    if fu>sf/2
        m=j-1;break
    end
    b=zeros(1,nfft);
    b(nl:nu)=a(nl:nu);
    b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
    c=ifft(b,nfft);
    yc(j)=sqrt(var(real(c(1:n))));
end
plot(fc(1:m),yc(1:m));


半个微笑 发表于 2015-11-11 17:05:44
楼主你好,你用的是什么声级计能导出WAV格式的音频

甫瑞德姆 发表于 2015-11-12 19:33:47
为什么计算声压有效值用 yc(j)=sqrt(var(real(c(1:nnn))));?

纪律先生 发表于 2015-11-29 21:59:55
这样在频域截取一段数据直接IFFT的方法简单粗暴,会引起时域数据波动,不可取。

为何不在频域直接积分得到结果?

纪律先生 发表于 2015-12-3 22:52:18
试了一下,代码中的IFFT有问题,结果的虚部很大,显然不正常。


simplebinbin 发表于 2015-12-21 16:11:31
半个微笑 发表于 2015-11-11 17:05
楼主你好,你用的是什么声级计能导出WAV格式的音频

bk声级计

simplebinbin 发表于 2015-12-21 16:12:35
甫瑞德姆 发表于 2015-11-12 19:33
为什么计算声压有效值用 yc(j)=sqrt(var(real(c(1:nnn))));?

ifft分析后含有虚部,只实部进行计算

simplebinbin 发表于 2015-12-21 16:13:34
纪律先生 发表于 2015-12-3 22:52
试了一下,代码中的IFFT有问题,结果的虚部很大,显然不正常。

先生可以把改进程序贴出来与大家一起分享

simplebinbin 发表于 2015-12-21 16:27:19
hebutkappa 发表于 2015-11-9 16:14
请问宋老师:
我有一个频域-声压曲线。横坐标HZ纵坐标db
可否直接用下面程序转成1/3倍频啊~做出来的数值 ...

横坐标频率分辨率怎么样,分辨率高的话,直接各个1/3倍频程频带内的db值计算有效值,分辨率低的话,先插值,在计算有效值
如有更好的方法,大家可以贴出来

大海1992 发表于 2016-1-13 14:00:09
学习了,真心不能闭门造车啊~~

JoofeeElec 发表于 2016-1-28 17:23:48
有点看不太懂,有数学模型吗?

simplebinbin 发表于 2016-1-29 13:35:10
JoofeeElec 发表于 2016-1-28 17:23
有点看不太懂,有数学模型吗?

程序已经很清楚了哦

songzy41 发表于 2016-1-29 21:26:29
1,有关LZ用了13633进行能量修正,可以参看:
https://www.ilovematlab.cn/thread-35655-1-1.html
2,IFFT变换后虚部数值很大,所以质疑IFFT有问题(39层),确实IFFT的设置中有问题。LZ和35层的程序都是基本上按王济的一书给出的,但该程序在IFFT中有问题。在负频率方面程序中为:
b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
应改为
b(nfft-nu+2:nfft-nl+2)=a(nfft-nu+2:nfft-nl+2);
3,为什么用 yc(j)=sqrt(var(real(c(1:n))));
IFFT后应该得到为实数,但是由于计算机有限字长的影响可能会有虚部,但数值很小。所以IIFT后取实部。在计算某频带的声压时应按均方根值来计算,var( .)是计算方差值,但因在某一频带内均方值和计算方差值相同(而均方值没有现成的函数),所以用var( .)来替代。再开根就得均方根值。

jwsj 发表于 2016-2-25 02:39:45

我想请问楼主,
您在15楼帖子中的第27行程序:xx=1.633*x.*w;         %加汉宁窗(恢复系数为1.633)

好像不对吧?运行会出错。原因是x是一个1行*n列的数组,而w是一个n行*1列的数组。如果你用数组乘'' .* '',数组维度不一致程序就会报错。

所以请问这里是该用矩阵乘:“xx=1.633*x*w ” ?  ; 或是该用w的转置进行运算 “xx=1.633*x.*(w.')”?;
亦或是其他方法?
我把以上两种猜测都试过了,结果相差很大。

simplebinbin 发表于 2016-5-9 15:45:25
jwsj 发表于 2016-2-25 02:39
我想请问楼主,
您在15楼帖子中的第27行程序:xx=1.633*x.*w;         %加汉宁窗(恢复系数为1.633)

使用xx=1.633*x.*(w')”;

孤独夜徘徊 发表于 2016-5-11 10:52:36
学习了,谢谢!
您需要登录后才可以回帖 登录 | 注册

本版积分规则

热门教程
站长推荐
快速回复 返回顶部 返回列表