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

  [复制链接]
mark8906 发表于 2018-1-12 10:15:29
本帖最后由 mark8906 于 2018-1-12 16:30 编辑
songzy41 发表于 2018-1-10 20:32
1,对于脉冲声可以对时间取平均
2,可以加时间计权

嗯。
宋老师,您看下面这两种情况我的理解对吗?
1,     48K采样频率,采样4096个点计算F时间计权的A声级。时间常数125ms,基本时间0.25ms,权重系数a = 0.002. 正常情况是每12个时域信号(频率计权后的信号)做为一个片段对每个时域信号的平方加时间计权求和算出这个片段的声压平方和。但是如果按这个帖子提供的计算声压级的方法,时域信号并没有提前做频率计权而是求1/3倍频声压级后进行了计权修正。
时间计权是对瞬时声压平方和的加权,而1/3倍频程需要先对信号进行FFT变换后取相应频带信号添0再FFT逆变换到时域再计算这个频带总的有效声压。那时间计权是对FFT反变换之后的信号进行平方加权求和,而不需计算有效声压?那计算量不是更大吗(有个论文说计算复杂度变小了,虽然循环次数变多了,不是很理解这个意思),并且这也不是真实的瞬时声压时域信号吧,可以这样吗?

2,   对噪声信号统一按48K采样的话,频率分辨率是不是很大了,还能检出较低频的中心频率1/3频带吗?还是说单纯计算声压级的话可以不需考虑这个,如果分析频谱功率分布的话再进行降采样处理?
3,   截取的一段频谱数据加多少长度的0算有意义,如果频谱只有1个长度有效数据,加0到4096逆变换回去这个还有意义没?

4,   这个帖子中的程序为什么不加低通滤波器,然后对原始信号先滤波处理呢?

这是2005年的一篇论文实现的时间计权,它printf函数里为什么将计算的声级减去7.8呢? ...

这是2005年的一篇论文实现的时间计权,它printf函数里为什么将计算的声级减去7.8呢? ...

zhr1111 发表于 2018-1-17 16:13:45
songzy41 发表于 2017-12-28 21:09
既然采集了1秒样本,对于噪声来说用统计的方法更好,所以建议用 pwelch函数计算功率谱,其中把nfft设为40 ...

宋老师  这里的f、f1、fc、cf是如何确定的啊?有没有具体的公式  看了这么多讨论  其实就没人能说清这个问题,而且为什么fc的第一个值是f1,麻烦您给解答一下
QQ图片20180117160908.png

zhr1111 发表于 2018-1-17 16:15:47
simplebinbin 发表于 2017-5-15 08:51
BK导出的数据也是进行处理了的,前面的帖子有说明。对于你的归一化后求传函这个问题,不是很懂。正常求没 ...

这里的f、f1、fc、cf是如何确定的啊?有没有具体的公式  看了这么多讨论  其实就没人能说清这个问题,而且为什么fc的第一个值是f1,麻烦您给解答一下,程序里这段是参考什么写的
QQ图片20180117160908.png

songzy41 发表于 2018-1-17 19:56:31
zhr1111 发表于 2018-1-17 16:15
这里的f、f1、fc、cf是如何确定的啊?有没有具体的公式  看了这么多讨论  其实就没人能说清这个问题,而 ...

在第2层给出的f、f1、fc、cf已作了注释:
f=[1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0]; %基准中心频率
f1=[20.00 25.0 31.5 40.0 50.0 63.0 80];
fc=[f1,100*f,1000*f,10000*f]; %%%%%%%%%中心频率%%%%%%%%
%20-16000Hz A声级计权值
cf=[-50.5,-44.7,-39.4,-34.6,-30.2,-26.2,-22.5,-19.1,-16.1,-13.4,-10.9,-8.6,-6.6,-4.8,-3.2,-1.9,-0.8,0,0.6,1.0,1.2,1.3,1.2,1.0,0.5,-0.1,-1.1,-2.5,-4.3,-6.6];
f和f1是基准中心频率,fc是中心频率,cf是A声级计权值。
1/3滤波器的中心频率是有国标的,可以查一下声学手册或实接查国标;A声级计权值也是有国标的。

zhr1111 发表于 2018-1-17 21:41:48
songzy41 发表于 2018-1-17 19:56
在第2层给出的f、f1、fc、cf已作了注释:
f=[1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0]; %基 ...

宋老师,想再问一下,这个国标的国标号具体是多少?我现在想做1/6倍频程,如果有相关频率规定我就不自己算了,因为算出来的值跟给定的貌似有差别。非常感谢!

songzy41 发表于 2018-1-18 10:25:31
本帖最后由 songzy41 于 2018-1-18 14:37 编辑
zhr1111 发表于 2018-1-17 21:41
宋老师,想再问一下,这个国标的国标号具体是多少?我现在想做1/6倍频程,如果有相关频率规定我就不自己 ...

我退休近20年了,许多资料已记不得了。我在文献中查到:
中华人民共和国国家质量监督检验检疫总局,中国国家标准化管理委员会.JJG 449-2001倍频程和1/3倍频程滤波器[S].北京:中国标准出版社,2005
不妨去查一下。我查到GB标准:GB-3241-1998-倍频程和分数倍频程滤波器

zhr1111 发表于 2018-1-18 11:08:20
songzy41 发表于 2018-1-18 10:25
我退休近20年了,许多资料已记不得了。我在文献中查到:
中华人民共和国国家质量监督检验检疫总局,中国 ...

谢谢您~

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

如果不进行傅里叶逆变换的话,直接在频域计算该频段内的有效声压级可以吗?

z527883588 发表于 2018-1-20 11:17:35
本帖最后由 z527883588 于 2018-1-20 11:34 编辑
songzy41 发表于 2016-1-29 21:26
1,有关LZ用了13633进行能量修正,可以参看:
https://www.ilovematlab.cn/thread-35655-1-1.html
2,IFFT变 ...

宋老师你好,我用的几个商业软件对有效值的定义并不是按均方根值来计算的,如图片所示。求方差时需要减去平均值然后再平方。
微信截图_20180120111530.png

songzy41 发表于 2018-1-21 19:46:35
z527883588 发表于 2018-1-20 11:17
宋老师你好,我用的几个商业软件对有效值的定义并不是按均方根值来计算的,如图片所示。求方差时需要减去 ...

对,用提供的算式更合理。

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

楼主有没有试过直接在频域内求该中心频率的有效值?

zhr1111 发表于 2018-1-22 15:21:06
songzy41 发表于 2018-1-18 10:25
我退休近20年了,许多资料已记不得了。我在文献中查到:
中华人民共和国国家质量监督检验检疫总局,中国国 ...

宋老师,再请教个问题,原始音频采样率32000hz,对应频率2000hz~3000hz频率段的频谱图需要单独处理,现在做降采样处理把采样率降低到8000hz,原来要研究的2000hz~3000hz频率段是不是就变到了500hz~750hz频率段了?

songzy41 发表于 2018-1-22 21:56:17
zhr1111 发表于 2018-1-22 15:21
宋老师,再请教个问题,原始音频采样率32000hz,对应频率2000hz~3000hz频率段的频谱图需要单独处理,现在 ...

通过降采样后需要处理的频率还是在2000hz~3000hz频率段,因为降采样仅把4000hz以上的分量滤除了。

zhr1111 发表于 2018-1-22 22:13:38
songzy41 发表于 2018-1-22 21:56
通过降采样后需要处理的频率还是在2000hz~3000hz频率段,因为降采样仅把4000hz以上的分量滤除了。
...

谢谢您~

simplebinbin 发表于 2018-2-27 10:35:32
z527883588 发表于 2018-1-20 10:06
如果不进行傅里叶逆变换的话,直接在频域计算该频段内的有效声压级可以吗? ...

求声压均方根时,定义在时域内进行。频域内计算我没这样做过,不太了解其具体意义。商业软件内,我试了一下,FFT分析后有线性和db选项,来显示fft后的幅值,用的20log(p/p0)这个公式,说明在频域内也是可以用的,感觉为了幅值显示方便,具体意义应该和时域内不一样。

Beau 发表于 2018-3-4 21:27:42
刚接触这个软件,看到这个帖子收获很多,感谢各位!
想请教一个问题,现在data=xlsread('G:\降噪-减振\轮胎试验2.1\Test2.xlsx',1,'A11504:C15009');%第1列时间,第2列声压(pa)
y=data(:,2); %db
x=data(:,1);%截取需要处理的数据段
运行时:
索引超出矩阵维度。
出错 Untitled (line 64)
Lp1(i)=20*log10(yc(i)/p0);%未计权1/3倍频程声压级
使声压级计算出问题了吧?这个公式是LP=20lg(p/p0);  
还有yc(j)=sqrt(var(real(c(1:n)))); 这一句没看明白。

Beau 发表于 2018-3-6 09:23:58
本帖最后由 Beau 于 2018-3-6 09:28 编辑

将代码调整如下,出来图1。
%A计权声压级频谱分析
clc;
clear;
close all;
data=xlsread('E:\Beau\Test2.xlsx',1,'A11504:C15009');
y=data(:,2); %db

fs=2000;%采样频率
p0=2e-5;%参考声压

f=[1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0]; %基准中心频率
f1=[20.00 25.0 31.5 40.0 50.0 63.0 80];
fc=[f1,100*f,1000*f,10000*f]; %%%%%%%%%中心频率%%%%%%%%%
%20-16000Hz A声级计权值
cf=[-50.5,-44.7,-39.4,-34.6,-30.2,-26.2,-22.5,-19.1,-16.1,-13.4,-10.9,-8.6,-6.6,-4.8,-3.2,-1.9,-0.8,0,0.6,1.0,1.2,1.3,1.2,1.0,0.5,-0.1,-1.1,-2.5,-4.3,-6.6];
t1=1;
t2=2;

x=data(:,1);%截取需要处理的数据段
n=length(x);
t=(0:1/fs:(n-1)/fs);

subplot(221);
plot(t,x);%瞬时声压时程图

w=hanning(n); %汉宁窗
xx=1.633*x.*w;%加汉宁窗(恢复系数为1.633)


nfft=2^nextpow2(n);
%nextpow2(n)-取最接近的较大2次幂
a = fft(xx,nfft);
f = fs/2*linspace(0,1,nfft/2);

w=2*abs(a(1:nfft/2)/n);
subplot(222);
plot(f,w');%绘制频谱图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1/3倍频程计算
oc6=2^(1/6);
nc=length(cf);
%下面这个求1/3倍频程的程序是按照振动振级计算那个来的
for j=1:nc
fl=fc(j)/oc6;
fu=fc(j)*oc6;
nl=round(fl*nfft/fs+1);
nu=round(fu*nfft/fs+1);
if fu>fs/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

aj_sumn=0;
for i=1:m
Lp1(i)=20*log10(yc(i)/p0);%未计权1/3倍频程声压级
end

%%%%%
for j=1:m
aj_sumn=aj_sumn+10^(0.1*Lp1(j));
end
Lp=10*log10(aj_sumn);%未计权总声压级
subplot(223);%绘制未计权1/3倍频程声压级图谱

bar(Lp1(1:m));
gg=zeros(1,m);
for i=1:m
gg(1:m)=fc(1:m);
end
ggg=1:m;
set(gca,'xtick',ggg);
set(gca,'xticklabel',gg);

%%%%%A计权1/3倍频程声压级
Lap=Lp1+cf(1:m);
aj_sum=0;
for j=1:m
aj_sum=aj_sum+10^(0.1*Lap(j));
end
LA=10*log10(aj_sum);%Aa计权总声压级

subplot(224);%绘制A计权1/3倍频程声压级图谱
bar(Lap(1:m));
gg=zeros(1,m);
for i=1:m
gg(1:m)=fc(1:m);
end
ggg=1:m;
set(gca,'xtick',ggg);
set(gca,'xticklabel',gg);
请问,cf 和 fu>fs/2 两个计数不一致导致频率只到了800吗?还是因为采样频率低(2000)?


图1

图1

Test2.xlsx

757.63 KB, 下载次数: 6


simplebinbin 发表于 2018-3-9 11:51:19
Beau 发表于 2018-3-6 09:23
将代码调整如下,出来图1。
%A计权声压级频谱分析
clc;

采样频率2000Hz,是能分析到800Hz的中心频率处,如要分析更宽的频带,请提高采样频率

Beau 发表于 2018-3-13 22:42:39
simplebinbin 发表于 2018-3-9 11:51
采样频率2000Hz,是能分析到800Hz的中心频率处,如要分析更宽的频带,请提高采样频率
...

好的,谢谢!:handshake

hhw527 发表于 2018-4-14 12:50:05
simplebinbin 发表于 2014-5-21 10:54
根据这两种方法计算结果还是有一些差别的,前一种方法称为方法1,后一种为方法2。计算abc.wav数据中前4s时 ...

楼主您好,那这个原始声压数据A计权以后应该是多少dB呢,0dB吗?

simplebinbin 发表于 2018-4-23 13:33:02
hhw527 发表于 2018-4-14 12:50
楼主您好,那这个原始声压数据A计权以后应该是多少dB呢,0dB吗?

当时用声级计随便采集了一段数据分析,原始数据结果多少忘了,这里主要讨论计算方法

hhw527 发表于 2018-4-24 16:34:59
simplebinbin 发表于 2015-1-13 10:47
最近在分析wav数据时发现,wav文件一般都做了归一化处理,MATLAB读取时需要进行修正(利用校准进行幅值修正 ...

楼主,我想问一下怎么修正,正常用这个原始声压数据,得到的A计权声级应该是0dB对吗

simplebinbin 发表于 2018-5-2 16:28:14
hhw527 发表于 2018-4-24 16:34
楼主,我想问一下怎么修正,正常用这个原始声压数据,得到的A计权声级应该是0dB对吗
...

1、原始数据对应A声级也不是0dB;
2、这个修正系数的技术需要一个标准的声校准器,用声级计采集一段标准声校准器发射的94dB标准声,导出时域数据,乘以一个修正系数,用这个程序刚好能计算到94dB,就是你需要的修正系数,不同的灵敏度,修正系数可能稍有偏差,每次测量计量都采集一段标准声对比。

一宝 发表于 2018-11-11 18:23:05
我想请问一下:程序中53、54行为什么那样把a赋值给b,直接对a进行ifft就可以了啊?就算幅值也可以直接写成b=a啊?

jiwei992112 发表于 2019-1-15 12:40:53
simplebinbin 发表于 2018-5-2 16:28
1、原始数据对应A声级也不是0dB;
2、这个修正系数的技术需要一个标准的声校准器,用声级计采集一段标准 ...

楼主您好,我是matlab小白,但是论文有个地方需要A计权1/3倍频程声压级图谱,我现在用软件计算模型得到了声压的时域关系数据(文件格式是.out,可以转换为txt格式),请问可以用你的那个第一个方法的程序计算么?  需要做些修改么? 我的采样频率是5000
您需要登录后才可以回帖 登录 | 注册

本版积分规则

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