plot fails with -Inf

9 次查看(过去 30 天)
Jonathan
Jonathan 2020-3-23
评论: Jonathan 2020-3-23
I need to plot a fft of an audio file. It appears that if any of the fft bins are zero, the conversion to dB produces an -Inf. Then the plot function does draw the figure but does not plot any data.
Is there an efficent way to trap the 0 value bins in the conversion to dB? For my work -Inf PSD would be equiv to zero.
The test audio file was a code generated audio file, 100Hz Sine for 10 seconds. Since only one fft bin will have a value, all the other bins are zero.
Thanks.
%Code Fragment:
[yS, Fs] = audioread(PathFN,'double');
L = length(yS);
Y = fft(yS,L);
% Version One
Px = Y.*conj(Y); % PSD Power of each freq components
Px = Px(1:L/2);
Px = 10*log10(Px); % Convert to dB
fVals = Fs * (0:(L/2)-1)/L;
plot(fVals, Px); % Plot
  3 个评论
Jonathan
Jonathan 2020-3-23
Hi, thank you for the response. The problem is that the plot function seems to failing because some of the y data is -inf
Peng Li
Peng Li 2020-3-23
did you try my solution? If you add a small number like eps to Px, there won't be any -inf.

请先登录,再进行评论。

回答(1 个)

Stijn Haenen
Stijn Haenen 2020-3-23
you can use Px(isinf(abs(Px)))=0; to set all -Inf values to 0
  3 个评论
Stijn Haenen
Stijn Haenen 2020-3-23
编辑:Stijn Haenen 2020-3-23
did you place it after the convertion to dB?
like this:
Px = Y.*conj(Y); % PSD Power of each freq components
Px = Px(1:L/2);
Px = 10*log10(Px); % Convert to dB
Px(isinf(abs(Px)))=0;
fVals = Fs * (0:(L/2)-1)/L;
plot(fVals, Px); % Plot
Does your Px also contains non real values?
Jonathan
Jonathan 2020-3-23
Yes, it was placed there and no just real values. aside from a typo, it does plot now. The single positive is at 100Hz as expected. However, it plots also at every 40Hz a spectra that is at/about -250db for the full range of Fs/2 except where the -Inf were. I didn't think it possible to have a negative PSD. Round off/precision errors from the program that generated the audio file? Thanks.

请先登录,再进行评论。

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by