필터 지우기
필터 지우기

Signal power with 1/3 octave filter

조회 수: 1 (최근 30일)
Luca Mastrangelo
Luca Mastrangelo 2017년 2월 9일
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).

답변 (0개)

카테고리

Help CenterFile Exchange에서 Measurements and Spatial Audio에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by