1/f amplitude spectrum and slope

Hi,

I would like to compute 1/f amplitude spectrum and loglog plot Or log10 .wav files and here is the code which does not work for many reasons. Any help is appreciated.

clear all;

close all;

clc;

[x,Fs] = audioread('*.wav');

Pxx = zeros(129,size(x,2));

for nn = 1:size(x,2)

[Pxx,F] = pwelch(x(:,nn),hamming(128),[],[],Fs,'psd');

end

PSD_x = 1./F(1:end).^2;

figure;

plot(log10(F),10*log10(PSD_x),'r','linewidth',4);

hold on;

plot(log10(F),10*log10(Pxx));

xlabel('loglog(Hz)'); ylabel('loglog(dB)');

Answer 1

[Mathieu NOE]

hello Mercede

Slightly modified code, simply give it a valid wav file name - it works as can be seen below - or where is your issue ?

clear all;
close all;
clc;
 
[x,Fs] = audioread('unfiltered_sound.wav'); % my test file name
 
% Pxx = zeros(129,size(x,2));
for nn = 1:size(x,2)
[Pxx,F] = pwelch(x(:,nn),hamming(128),[],[],Fs,'psd');
end
 
PSD_x = 1./F(1:end).^2;
figure;
% plot(log10(F),10*log10(PSD_x),'r','linewidth',4);
semilogx(F,10*log10(PSD_x),'r','linewidth',4);
hold on;
semilogx(F,10*log10(Pxx));
% plot(log10(F),10*log10(Pxx));
% xlabel('loglog(Hz)'); ylabel('loglog(dB)');
xlabel('(Hz)'); ylabel('(dB)');

[Mercede Erfanian]

Thank you for the quick response. Yes the plot works for me as well. However, it seems the angle of the slope does not change when I run it on different .wav files. I expect the slope to be very shallow for whitenoise as oppose to bird chirpping (with high frequency content). I would also like to print the alpha value of (1/f) of each sound. Would you please help?

Thank you very much

[Mathieu NOE]

Ok , so the 1/f noise slope is not truly a fixed slope curve (as I first understood) , does it mean you want to generate a loglog linear slope / average curve (linear in log log coordinates). Otherwise, your code as it is today will not change it slope according to different wav files.

that would mean a kind of fitting the psd data with a loglog linear curve ?

maybe I should try on a couple of your wav files (if you can share them ?)

[Mercede Erfanian]

-Ok , so the 1/f noise slope is not truly a fixed slope curve (as I first understood) = Exactly, the slope must change with each sound and its acoustic characteristics (please see the screenshot 1 attached).

-does it mean you want to generate a loglog linear slope / average curve (linear in log log coordinates). Otherwise, your code as it is today will not change it slope according to different wav files. = Yes, you are right, actually loglog plot with line of best fit is a better option than semilogx (please see the screenshot 2 attached).

-Would be also great to be able to print the value of alpha with the plot.

I have also attached a few .wav file.

Thank you very much for the help,

[Mercede Erfanian]

Here are the sounds

[Mercede Erfanian]

Basically what I need is pretty much like this example but I have hard time adapting it to run on a .wav file.

https://www.mathworks.com/matlabcentral/answers/473754-plot-best-fit-line-in-log-space-get-slope

[Mathieu NOE]

hello again

so this is a new code including the loglog fit submission above

I tested it on the four wav files - seems to do the job !

hope it helps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% data
[data,Fs] = audioread('2.wav');

channel = 1;
signal = data(:,channel);
samples = length(signal);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000;    % 
OVERLAP = 0.75;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);
% remove data for freq = 0
ind = find(freq>0);
freq = freq(ind);
sensor_spectrum = sensor_spectrum(ind);


%% log log poly fit

Bp = polyfit(log10(freq), log10(sensor_spectrum), 1);
Yp = polyval(Bp, log10(freq));
Yp2 = 10.^(Yp);
%

figure(1),loglog(freq,sensor_spectrum,'b',freq,Yp2,'r');grid on
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude');
expstr = @(x) [x(:).*10.^ceil(-log10(abs(x(:))))  floor(log10(abs(x(:))))]; % Mantissa & Exponent
Bp2 = expstr(10^Bp(2));
eqn = sprintf('y = %.2f\\times10^{%d}\\cdotx^{%.2f}',Bp2(1), Bp2(2), Bp(1));
legend('Data', eqn)
text(freq(2), sensor_spectrum(end) , strcat('Slope in Log-Log Space =', num2str(Bp(1))), 'Interpreter', 'none');

function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select);
freq_vector = (select - 1)*Fs/nfft;
end 

[Mercede Erfanian]

Dear Mathieu,

I can not thank you enough. I really appreciate your time and help,

Cheers.

[Mathieu NOE]

You're welcome !

[ANURAG SINGH]

you're good