A计权声压级的计算（matlab)

# 该文件tiaoshi.m记录了一个声卡测试声压的简易过程
# 主要完成了无计权和A计权两种状态的声压级
# 未来有时间可以增加的功能：1. 加窗fast（125ms),slow(1s) 计算 连续声压级
# 连续瞬时成分分析


# 用到的材料和数据：
# 1. 1Ksine.wav  用于检验计权对频率和幅度的影响。不管什么计权，1K音源的无计权和计权声压应该一致，
# 其频谱中显示的频率和幅度两种状态下应该一致。
# 2. 1k94dBRecording副本av 以及1k94dBRecording滤波.wav
# 前者是4231的录音，本身有一些环境噪声，后者是对前者的滤波处理，尽可能只保留1KHz左右的频率
#相关结论： a。环境噪声对声压级测试影响不大。b.这个音源的声压明确为94dB，无论无计权或者A计权
#3. whitelongrecording.wav   白噪声的录音 声压级85.4
#4.lonely.wav  录音 声压计测试值 82.9
#5. wale wale.wav   录音 声压计测试值 83.4
#6. califorlia.wav 录音 声压计测试值 80.8
#3~6 对比音频文件录音算出来的值和声压级测试值证明程序可行。


# 麦克风灵敏度 22.78mv/Pa
# 声卡 Vpeak = 1.8
# 前置放大器增益放到0dB，就是放大器对麦克风的输出电平不做处理


实际的一些问题：
a. 声卡标准input有效值4dBV（1.58V),最大input有效值11dBV(3.54)
    但是厂商设置的Vpeak =1.8V
    dBV 和电压转换 dBV = 20log（电压/1V)
b. 麦克灵敏度50mv/pa， 校准后22.78mv/pa.   跟上次校准（33mv/Pa）还是差异比较大。
c. 最早的尝试失败的原因：【麦克风灵敏度和声卡峰值电平与规格书不同，需要采用实际校准值；】
                         【先求了瞬时声压序列，然后对声压做FFT，最后发现频率成分全错了。
                         改为先对电平序列做fft，ifft处理，加权完成后再去求声压，是正确的，
                         问题应该出在复合函数的FFT变换上；】
                         【fft时对幅度做了归一化处理，然后ifft回去的时候就很尴尬，最后选择这两部都不做归一化处理】






d. 每一步应该看中间值，以及plot出来看频率和幅度。


处理顺序：
1. 读取录音文件
2. 求无计权声压：直接推到声压序列，求序列RMS值，得到无计权声压
2. A计权声压： 
    对录音文件的rawdata进行fft处理：
    乘以A计权网络
    ifft回去，得到加权过的rawdata
    推到声压序列
    求得A计权声压

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

# 各步骤解释


[x,fs] = audioread('1Krecordingwav');   %替换成其他录音


L = length(x);
y=fft(x);     # 注意：没有做归一化处理
f  =fs/2.*linspace(0,1,L/2+1);
y = y(1:L/2+1);
plot(f,abs(y),'b'); % 读取wav中的信息，并做FFT变换。 检查是否只有1K，以及1K的幅度(0.15)，最大频率。
% 到此时，只有1K，幅度为0.15， 最大频率2万五。


A= filterA(f);
yweighted = A'.*y; %得到加权后的complex Amplitude
plot(f,abs(yweighted),'r'); % 这个时候一切正常，仍然只有1K，1K的幅度为0.15，最大两万五


%然后准备ifft变换
xweighted = ifft(yweighted); %得到和yweighted 相同长度的xweighted # 注意：没有做归一化处理
xweighted = real(xweighted);  % 处理后的xweighted不在是一个实数序列，只取其实部
%检验一下xweighted的内容
tempdata = xweighted;
tempL = length(tempdata);
temp1 = fft(tempdata)/tempL*2;
plot(f,abs(temp1),'g');
% 这个时候也没有问题，注意的是temp1的运算，不再需要除以L再乘以2.因为你已经在实际幅度上运算了。
% 也许前面不需要做归一化，归一化应该放到最后面。%


# x， xweighted 分别代表了无计权的原始数据和A计权处理后的数据。
#下一步根据数据去推声压，
#这一部分是求A计权声压
linelevel =1.8;
preampout = abs(xweighted.*1.8);
preampin = preampout;
pressure = 1000.*preampin./22.78;
pref = 0.00002;
pressure = sqrt(mean(pressure.^2));
AweightedSPL  =20.*log10(pressure/pref);


#这一部分是求无计权声压
pshorttime = abs(x.*1.8);
pressureshorttime = 1000.*pshorttime./22.78;
pressureshorttime = sqrt(mean(pressureshorttime.^2));
SPL  =20.*log10(pressureshorttime/pref);


#返回无计权和A计权两种声压值
SPL
AweightedSPL