Increasing value over time of Power of a audio signal

Hi everyone, i have an audio signal and the power of the same signal increases over time (with more samples) so i have a very big value of the power and a wrong power spectrum if i use this code

% Lunghezza del segnale e trasformata FFT
m = length(xn);                 % lunghezza del segnale originale

Unrecognized function or variable 'xn'.

y = fft(xn, m);                 % calcolo della FFT del segnale

% Vettore delle frequenze
f = (0:m-1) * (fs / m);           % frequenze associate ai campioni FFT

% Calcolo dello spettro di potenza normalizzato
power = abs(y).^2 / m;            % spettro di potenza normalizzato

% Plot dello spettro di potenza
figure;
plot(f(1:floor(m/2)), power(1:floor(m/2)));  % plottiamo solo la metà positiva
xlabel('Frequency (Hz)');
ylabel('Power');
title('Power Spectrum of Signal');

But if i use the function pspectrum or the signal analyzer app i have different results much smaller

[p, f] = pspectrum(xn, fsn, 'power'); % Calcola la potenza spettrale e la frequenza
plot(f, p); % Plotta la potenza spettrale in funzione della frequenza
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');
title('Power Spectrum of Signal');

same plot for pspectrum and signal analyzer so i think it s correct, why i have that error with power over 160 db/hz? i need to calculate the power myself cause i need the values for a model thank you

Answer 1

[William Rose]

@Alex Kouroupis,

You say "I need to calculate the power myself cause i need the values for a model".

The first plot you provided has the vertical axis label cut off, so I am not sure how it differs from the second plot. Normalization of the power spectrum depends on several choices, including whether you want the power spectrum or the power spectral density, one-sided versus two-sided, winodw choice and possible adjustment for windowing, etc. See here for a discussion of these choices, and how to get results from fft() that match the results from other matlab approaches to estimating the spectrum. The second plot you show uses pspectrum() to compute the spectrum. pspectrum() makes a lot of decisions based on the length of the time-domain sequence, etc. See here for details. It will be a lot of work to write your own code that gives the same output as pspectrum.

If you want to compute the power yourself, then do it by a consistent and clear method, the same way every time. The values you obtain will work with the model you develop.

[William Rose]

@Alex Kouroupis,

[Matlab answers not letting me run code below. ("You sare not authorized to perform this operation.") It worked earlier. Maybe it will work for you. I hope there are no errors since I had to copy/paste some.]

Make three similar signals. All three have two sinusoids with same frequencies and amplitudes. They differ in sampling rate and duration

N1=50;  fs1=50;  t1=(0:N1-1)/fs1;  % constants
N2=100; fs2=100; t2=(0:N2-1)/fs2;  % constants
N3=50;  fs3=100; t3=(0:N3-1)/fs3;  % constants
fa=8; fb=16;  % frequencies of sinusoids
x1=cos(2*pi*fa*t1)+0.5*sin(2*pi*fb*t1); % time domain signal 1
x2=cos(2*pi*fa*t2)+0.5*sin(2*pi*fb*t2); % time domain signal 2
x3=cos(2*pi*fa*t3)+0.5*sin(2*pi*fb*t3); % time domain signal 3

Plot time domain signals

figure; subplot(211); plot(t1,x1,'-rx',t2,x2,'-go',t3,x3,'-b+')
grid on; xlabel('Time (s)'); ylabel('Ampltude'); 
legend('x1','x2','x3')

Compute power spectra with pspectrum.

[P1,f1]=pspectrum(x1,fs1,'power');
[P2,f2]=pspectrum(x2,fs2,'power');
[P3,f3]=pspectrum(x3,fs3,'power');

Display frequency spacing of the spectra:

fprintf('Delta f for pspectrum spectra 1,2,3 = %.4f, %.4f, %.4f Hz.\n',...
    f1(2)-f1(1),f2(2)-f1(1),f3(2)-f3(1));

Delta f for pspectrum spectra 1,2,3 = 0.0061, 0.0122, 0.0122 Hz.

Plot spectra

subplot(212); plot(f1,P1,'-rx',f2,P2,'-go',f3,P3,'-b+');
grid on; xlabel('Frequency (Hz)'); ylabel('Power'); 
legend('P1','P2','P3'); title('pspectrum')

Example above shows that pspectrum() has smoothed out the spectrum with a Kaiser window - otherwise the spectrum would have looked like two delta functions at 8 Hz and 16 Hz. It shows that the power returned by pspectrum() is not in dB. It shows that the power returned by pspectrum is not the power spectral density (psd). (If it were the psd, then the values at the peaks would look different by a factor of two when delta f changes by a factor of two.) Therefore you will want to adjust the y-axis labels on the plots in your original post.

More later.

[William Rose]

@Alex Kouroupis,

Related to previous comment:

N1=50;  fs1=50;  t1=(0:N1-1)/fs1;  % constants
N2=100; fs2=100; t2=(0:N2-1)/fs2;  % constants
N3=50;  fs3=100; t3=(0:N3-1)/fs3;  % constants
fa=8; fb=16;  % frequencies of sinusoids
x1=cos(2*pi*fa*t1)+0.5*sin(2*pi*fb*t1); % time domain signal 1
x2=cos(2*pi*fa*t2)+0.5*sin(2*pi*fb*t2); % time domain signal 2
x3=cos(2*pi*fa*t3)+0.5*sin(2*pi*fb*t3); % time domain signal 3

Plot time domain signals

figure; subplot(211); plot(t1,x1,'-rx',t2,x2,'-go',t3,x3,'-b+')
grid on; xlabel('Time (s)'); ylabel('Ampltude'); 
legend('x1','x2','x3')

Compute power spectra with fft().

P1=abs(fft(x1)/N1).^2; 
P1=[P1(1),P1(2:N1/2)*2,P1(N1/2+1)]; % convert 2 sided to 1-sided
f1=fs1*(0:N1/2)/N1;
P2=abs(fft(x2)/N2).^2; 
P2=[P2(1),P2(2:N2/2)*2,P2(N2/2+1)]; % convert 2 sided to 1-sided
f2=fs2*(0:N2/2)/N2;
P3=abs(fft(x3)/N3).^2; 
P3=[P3(1),P3(2:N3/2)*2,P3(N3/2+1)]; % convert 2 sided to 1-sided
f3=fs3*(0:N3/2)/N3;

Display frequency spacing of the spectra:

fprintf('Delta f for FFT spectra 1,2,3 = %.2f,%.2f,%.2f Hz.\n',...
    f1(2)-f1(1),f2(2)-f1(1),f3(2)-f3(1));

Delta f for FFT spectra 1,2,3 = 1.00,1.00,2.00 Hz.

Plot spectra

subplot(212); 
plot(f1,P1,'-rx',f2,P2,'-go',f3,P3,'-b+');
grid on; xlabel('Frequency (Hz)'); ylabel('Power'); 
legend('P1','P2','P3'); title('fft')

Example above shows that the powers computed with fft(), with adjustment as shown in the formulas I used, match the values returned by pspectrum() in my previous comment.

[William Rose]

@Alex Kouroupis,

The plot in my previous comment, and the plot in the comment before that, show power=0.5 at 8 Hz and power=0.125 at 16 Hz. These values exactly equal the rms power of the 8 Hz and 16 Hz sinusoids. If the sinusoids were voltages, across a 1 ohm resistor, that's the power in watts that they would produce.

Regard the spectra returned by pspectrum() with a critical eye. The frequencies are spaced at 0.006 or 0.012 Hz. I think this is misleading, because the signals are only 0.5 or 1 second long. With signals of such duration, it is not possible to distinguish frequencies so close together. Also, the sum of the spectrum values, with pspectrum(), does not equal the total power in the signal. Whereas for the spectrum we got with abs(fft(x)/N).^2, the sum of the spectrum values does equal the total power of the signal. On the other hand, if the signal did not happen to be an exact integer number of cyles during the period sampled, the spectrum from pspectrum would be superior in terms of leakage to the spectrum from fft(), since pspectrum() uses a window to control leakage.

[Alex Kouroupis]

@William Rose

Thank you very much for the answer apreccieated. So is it safe to get the power of a signal with pspectrum() function? I need to do some clssifications of a signal and i use power and energy as featues (with others), if i put more samples in as frequencyresolution will i get a better pspectrum to see near frequencies? also how many samples does pspectrum() use? can t find the default prametres of the fucntion. Thank you veryyyyy very much

[William Rose]

@Alex Kouroupis,

If you want signal power for classificaiton purposes, why not just compute the variance of the time domain signal?

I suspect that if you are consistent in your use of pspectrum() with various signals, then the outpout from pspectrum() will work just as well for classificaiton as if you used abs(fft()/N).^2 or cpsd() or some other approach.

pspectrum() uses as many samples as you pass to it. It appears, from the spectra produced by pspectrum(x,fs,'power'), that it pads with zeros. I agree that it is not obvious from the documentation.