0

I have filtered my wav sound in a specific band with the butterworth bandpass filter. Now i want to calculate the intensity of that sound in the specific band.

[B1,A1] = butter(4,[10 635]/(fs_orig/2),'bandpass');
freqz(B1,A1) %The frequency response of your filter1
dataIn = x1; %your music
dataOut1 = filter(B1, A1, dataIn); % // filter command filters your music

So first have to calculate the rms value of the specific band? rms1 = rms(dataOut1) this gives me rms between 10 and 625 hz. But how should I calculate it into dB when 0 db corresponding to an rms value of 1?

mehmet
  • 1,631
  • 16
  • 21
user4390280
  • 33
  • 10
  • If i want to plot that into the frequency domain i should use plot(abs(fft(20*log10(rms1))));xlabel('Frequency(Hz)');ylabel('Amplitude'); is that correct? – user4390280 Jan 06 '15 at 20:21
  • You normally take fft of amplitude or power and take log after doing fft. – Navan Jan 06 '15 at 20:29
  • You mean like this? rmslog = 20*log10(rms1); fftvalue = fft(rmslog); %plot plot(abs(fftvalue);xlabel('Frequency(Hz)');ylabel('Amplitude'); – user4390280 Jan 06 '15 at 20:43
  • When i try to plot it doenst look like it has to be.. Dont know what im missing – user4390280 Jan 06 '15 at 22:06

1 Answers1

0

Your approach is good (including the 20*log10(rms), as suggested), but you will end up plotting the instantaneous amplitude...ie, the amplitude for each sample of your data, which is likely at a sample rate of 44100 Hz or something. This is probably not want you want because, as humans, this is not how we perceive "volume" or "loudness".

A human cannot hear amplitude changes that are faster than, say, 10 to 40 Hz. So, you need to lowpass filter your instantaneous amplitude values to something slower than your native 44100 Hz sample rate.

I recommend something like this:

%filter your data
[B1,A1] = butter(4,[10 635]/(fs_orig/2),'bandpass');
freqz(B1,A1) %The frequency response of your filter1
dataIn = x1; %your music
dataOut1 = filter(B1, A1, dataIn); % // filter command filters your music

%compute the instantaneous power
data_pow = dataOut1.^2;

%low pass filter at 5 Hz
lp_Hz = 5;
[b,a]=butter(2,lp_Hz/(fs_orig/2));
data_pow = filter(b,a,data_pow);

%convert to dB
data_dB = 10*log10(data_pow);

%plot
t_sec = ([1:length(data_dB)]-1)/fs_orig;
plot(t_sec,data_dB);
xlabel('Time (sec)');
ylabel('Amplitude (dB)');
chipaudette
  • 1,655
  • 10
  • 13
  • Thanks, but if i should take a better range like [1000 1625], how do i need to plot this amplitudes in db? – user4390280 Jan 07 '15 at 01:41
  • Assuming that you are talking about a frequency range of [1000 1625], you simply change the first line to be: `[B1, A1]=butter(4,[1000 1625]/(fs_orig/2),'bandpass');`. The rest of the code will be the same. The code already plots in dB. – chipaudette Jan 07 '15 at 12:09
  • How can i plot the frequency spectrum? – user4390280 Jan 07 '15 at 12:21
  • Whose spectrum? Your `freqz` command should have plotted the frequency spectrum of your filter. Alternatively, if you want the frequency spectrum of your original sound, that's a whole different question than what you asked. Search for "spectrum" or "fft" and you'll find lots of hits. If you've gotten the answer to your original question about plotting intensity in dB, you should consider marking this question "answered". – chipaudette Jan 07 '15 at 14:03
  • But when i plot it my power has a negative dB value, is that possible? – user4390280 Jan 08 '15 at 20:24
  • Also when i plot my sound with this code it shows me a negative db value and then it increases to a positive db in time. – user4390280 Jan 08 '15 at 20:32
  • Don't forget that "dB" is not the same as "dB SPL". Negative values perfectly normal for analyzing audio data, especially if you used a function like `wavread` to read them in. zero dB means that the rms = 1. A negative dB means that the rms is less than one. No problem. If, however, you want "dB SPL" (sound pressure level) you need a calibrated microphone. – chipaudette Jan 09 '15 at 21:56