计算声音OverAll图及1/3 倍频程图Mathematica编程...
本帖最后由 rdh 于 2019-5-26 01:05 编辑程序编了很长一段时间，现分享出来。（备注：程序不够精简，还有一些变量命名不统一，请大家将就一些看吧，提一点意见；现在正在编写直流电机，由电流信号计算输出转速信号的程序，有兴趣的可以一些探讨）此贴转载或复制，请注明出处或者作者rdh,谢谢！
先看程序准确度和效果。图1和图3为海德声科Head Acoustics 软件计算；图2和图4为Mathematica计算软件编程的结果；
程序如下：
fs = 48000;(*采样频率/Hz;*)
dt = 1/fs;(*每隔\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"][][[
15 ;; (ntotal + 14), 2]];
timedata1 =
TimeSeries}];
ListPlot
ts = dt n;
(*一个样本周期为T=dt n*)
loop = IntegerPart;(*总的帧数，或者数据块的个数*)
tlp = Table;
j1 = 0;
Do[
j1 = j1 + 1;
Print["第", j1, "帧样本 时间点t=", N[];
onesample =
timedata[]];(*FFt分析的一个样本;*)
onesample1 = TimeSeries}];
Print];
w = Table[
N], {x, -1/2, 1/2, 1/(
n - 1)}];(*HannWindow窗数据表形式恢复系数为1.633*)
newonesample = w onesample;
newonesample1 =
TimeSeries}];
Print];
(*傅里叶变换*)
a = Fourier;
fft = TimeSeries, {fs/n Range}];
Print[ListLinePlot[fft,
PlotRange -> {{0, fs/2}, {0, 0.15}},
PlotLabel -> "样本-频域信号"]];
(*%1/3倍频程计算*)
oc6 = 2^(1/6);
nc = Length;
yc = Table;
lp1 = Table;
j = 0;
Do[
j = j + 1;
fl = fc[]/oc6;
fu = fc[]*oc6;
nl = Round + 1];
nu = Round + 1];
b = Table;
b[]] = a[]];
b[]] =
a[]];
c = InverseFourier;
yc[] = Sqrt^2]];
lp1[] = 20*Log10]/p0];
, {l, 1, 30}];
lpa = 10 Log10];
(*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[] = 10 Log10];
Print["A计权后总的声压级：LP(A)", "=", tlp[], "dB"];
Print[BarChart[{lp1 + cf}, PlotLabel -> "A计权后1/3倍频程图",
LabelingFunction -> Above, ChartLabels -> axial]]
, {b1, 1, loop}];
overall = TimeSeries}];
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 发表于 2019-5-26 08:59
支持楼主的原创！
谢谢！ luckyhuman1 发表于 2019-5-26 13:12
这么厉害，膜拜
谢谢 这么厉害，膜拜 用Matlab更容易实现吧，我偶尔用用MMA，也是边学边用，用了MAT回来用MMA感觉好不舒服，不过系统学习MMA也是很有帮助的，技多不压身。 mxlzhenzhu 发表于 2019-6-2 20:54
用Matlab更容易实现吧，我偶尔用用MMA，也是边学边用，用了MAT回来用MMA感觉好不舒服，不过系统学习MMA也是 ...
个人感觉MMA很强大。
页:
[1]