计算声音OverAll图及1/3 倍频程图Mathematica编程...

rdh

程序编了很长一段时间，现分享出来。（备注：程序不够精简，还有一些变量命名不统一，请大家将就一些看吧，提一点意见；现在正在编写直流电机，由电流信号计算输出转速信号的程序，有兴趣的可以一些探讨）此贴转载或复制，请注明出处或者作者rdh,谢谢！
先看程序准确度和效果。

图1和图3为海德声科Head Acoustics 软件计算；图2和图4为Mathematica计算软件编程的结果；

程序如下：
fs = 48000;(*采样频率/Hz;*)
dt = 1/fs;(*每隔\[CapitalDelta]t，采集一个数据点*)
t = 10;(*采样时间/s;*)
ntotal = t/dt;(*总共采集的数据点个数;*)
n = 8192;(*spectrumsize/频谱数;*)
df =fs/n;(*频率分辨率/Hz;*)
p0 = 2*10^-5;(*参考声压/pa;*)
fc = {20.00, 25.0, 31.5, 40.0, 50.0, 63.0, 80, 100.`, 125.`, 160.`,
  200.`, 250.`, 315.`, 400.`, 500.`, 630.`, 800, 1000.`, 1250.`,
  1600.`, 2000.`, 2500.`, 3150.`, 4000.`, 5000.`, 6300.`, 8000, 10000,
   12500, 16000};
(*%1/3倍频程中心频率*)
tt = 50 %;(*重叠率为50%*)
(*time Single;*)
timedata =
Import["D:\\database\\##.xlsx"][[1]][[
  15 ;; (ntotal + 14), 2]];
timedata1 =
  TimeSeries[timedata, {dt Range[0, (ntotal - 1)]}];
ListPlot[timedata1, PlotLabel -> "采集的原始信号"]
ts = dt n;
(*一个样本周期为T=dt n*)
loop = IntegerPart[2 t/ts - 1];(*总的帧数，或者数据块的个数*)
tlp = Table[0, {g, loop}];
j1 = 0;
Do[
  j1 = j1 + 1;
  Print["第", j1, "帧样本   时间点t=", N[[dt (1 + j1  n/2)]];
  onesample =
   timedata[[Range[(1 + (j1 - 1) n/2), (j1 + 1) n/2]]];(*FFt分析的一个样本;*)


  onesample1 = TimeSeries[onesample, {dt Range[1, n]}];
  Print[ListPlot[onesample1, PlotLabel -> "截取分析的一个样本-时域信号"]];
  w = Table[
    N[1.633 HannWindow[x]], {x, -1/2, 1/2, 1/(
     n - 1)}];(*HannWindow窗数据表形式恢复系数为1.633*)
  newonesample = w onesample;
  newonesample1 =
   TimeSeries[newonesample, {dt Range[1, n]}];
  Print[ListPlot[newonesample1, PlotLabel -> "加汉宁窗处理后的样本-时域信号"]];
  (*傅里叶变换*)
  a = Fourier[newonesample];
  fft = TimeSeries[Abs[a], {fs/n Range[1, n]}];
  Print[ListLinePlot[fft,
    PlotRange -> {{0, fs/2}, {0, 0.15}},
    PlotLabel -> "样本-频域信号"]];
  (*%1/3倍频程计算*)
  oc6 = 2^(1/6);
  nc = Length[fc];
  yc = Table[0, {g, nc}];
  lp1 = Table[0, {g, nc}];
  j = 0;
  Do[
   j = j + 1;
   fl = fc[[j]]/oc6;
   fu = fc[[j]]*oc6;
   nl = Round[fl*n/Subscript[f, s] + 1];
   nu = Round[fu*n/Subscript[f, s] + 1];
   b = Table[0, {g, n}];
   b[[Range[nl, nu]]] = a[[Range[nl, nu]]];
   b[[Range[n - nu + 1, n - nl + 1]]] =
    a[[Range[n - nu + 1, n - nl + 1]]];
   c = InverseFourier;
   yc[[j]] = Sqrt[1/n Total[Re[c]^2]];
   lp1[[j]] = 20*Log10[yc[[j]]/p0];
   , {l, 1, 30}];
  lpa = 10 Log10[Total[10^(0.1 lp1)]];
  (*Print["总的声压级：LP","=",lpa,"dB"];*)
  axial = {"20.00", "25.0", "31.5", "40.0", "50.0", "63.0", "80",
    "100.`", "125.`", "160.`", "200.`", "250.`", "315.`", "400.`",
    "500.`", "630.`", "800", "1000.`", "1250.`", "1600.`", "2000.`",
    "2500.`", "3150.`", "4000.`", "5000.`", "6300.`", "8000", "10000",
     "12500", "16000"};
  Print[BarChart[lp1, PlotLabel -> "1/3倍频程图",
    LabelingFunction -> Above, ChartLabels -> axial]];
  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};
  tlp[[j1]] = 10 Log10[Total[10^(0.1 (lp1 + cf))]];
  Print["A计权后总的声压级：LP(A)", "=", tlp[[j1]], "dB"];
  Print[BarChart[{lp1 + cf}, PlotLabel -> "A计权后1/3倍频程图",
    LabelingFunction -> Above, ChartLabels -> axial]]
  , {b1, 1, loop}];
overall = TimeSeries[tlp, {ts/2  Range[1, loop]}];
Print[ListLinePlot[overall, PlotRange -> {{0, 10}, {15, 60}},
   PlotLabel -> "声压级的OverAll图像", PlotTheme -> "Detailed"]];


计算的过程详细分解如下：
第一步：采集的10秒的原始信号：

第二步：取第1帧数据块（有谱线数来决定）；第1帧样本   时间点t=0.0853542

第三步：对第1帧数据块加汉宁窗；

第四步：对第1帧数据块FFT变换；

第五步：由FFT变换的频率信号计算转化为1/3倍频图；

第六步：1/3倍频进行A计权，同时计算出该数据块的声压级；A计权后总的声压级：LP(A)=36.332dB

第七步：取第2帧数据块（由重叠率来决定），重复第二步带第七步；
一直取到把10秒的数据全分析完。（就是把10秒的数据分解成一份一份的数据块）

第八步：把各个数据块声压级按照时间顺序连成一条线则得到OverAll图；

好了，以上工作已经完成。

对于稳定的信号，测试10秒，最终出来一个声压级值和1/3倍频程图。它是由所有数据块的声压级值，线性平均后得到。

jdx_2007

支持楼主的原创！

luckyhuman1

这么厉害，膜拜

rdh

jdx_2007 发表于 2019-5-26 08:59
支持楼主的原创！

谢谢！

rdh

luckyhuman1 发表于 2019-5-26 13:12
这么厉害，膜拜

谢谢

wolfxin

这么厉害，膜拜

mxlzhenzhu

用Matlab更容易实现吧，我偶尔用用MMA，也是边学边用，用了MAT回来用MMA感觉好不舒服，不过系统学习MMA也是很有帮助的，技多不压身。

rdh

mxlzhenzhu 发表于 2019-6-2 20:54
用Matlab更容易实现吧，我偶尔用用MMA，也是边学边用，用了MAT回来用MMA感觉好不舒服，不过系统学习MMA也是 ...

个人感觉MMA很强大。