Effacer les filtres
Effacer les filtres

Signal power with 1/3 octave filter

1 vue (au cours des 30 derniers jours)
Luca Mastrangelo
Luca Mastrangelo le 9 Fév 2017
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).

Réponses (0)

Catégories

En savoir plus sur Measurements and Spatial Audio dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by