Signal power with 1/3 octave filter
2 次查看(过去 30 天)
显示 更早的评论
Hi, I'm trying to calculate the power in dB of a signal and express it in 1/3 octave. I mean, I apply the 1/3 octave filter to the signal and then calculate the power for each band. The problem is that when I compare my program with a professional one I get very different results. Below my code:
[x, Fs] = audioread('audio/pink/signal.wav');
C = 25;
%%octave
lx = length(x);
% filter design
Fc = [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 ...
500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 ...
6300 8000 10000 12500 16000 20000];% frequenze centrali ANSI
Fc = Fc(Fc<Fs/2);
fl = Fc*2^(-1/6); % lower freq
fu = Fc*2^(1/6); % upper freq
fu(fu>Fs/2) = Fs/2 -1;
numBands = length(Fc);
ord = 12;
z = zeros(2*ord, numBands);
p = zeros(2*ord, numBands);
k = zeros(1, numBands);
for n = 1:numBands
[z(:,n),p(:,n),k(n)] = cheby1(ord,10,[fl(n) fu(n)]/(Fs/2));
end
% apply filter
x_oct = zeros(lx, numBands);
s = zeros(ord, 6, numBands);
for i = 1:numBands
s(:, :, i) = zp2sos(z(:,i),p(:,i),k(i));
x_oct(:, i) = sosfilt(s(:,:,i),x);
end
x_pow_oct = zeros(1, numBands);
for k = 1:numBands
fftx = fft(x_oct(:,k));
fftxsquared = fftx'*fftx;
s1 = sum(fftxsquared/lx);
x_pow_oct(k) = 10*log10(s1) + C;
end
%%end function
maxV = max(x_pow_oct);
f = [0 Fc(5:5:30)];
figure(150)
subplot(1,2,1)
bar(x_pow_oct);
title('x_oct_pow');
set(gca,'XTickLabel',{round(f)});
xlabel('Central frequencies');
ylabel('Power');
ylim([10 maxV+10]);
xlim([0 numBands+1]);
To show the differences with the professional software, below my result (first image) and the professional one (seond image).
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/160531/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/160536/image.bmp)
0 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Measurements and Spatial Audio 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!