# 该文件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
