Extract the relative amplitude from various signals to a matrix for a powerspectrum plot

Hi! I am still new to MAtTLAB and hope that you can help med with something I'm struggeling with.

I want to plot a power spectrum (See below) so i can investigate the relative powerdistributions across frequencies for multiple signals that I am interested in. I have a list of start and endtimes of multiple sections from multiple wav-files which i have loaded into Matlab through this loop (dataX is the list):

selpath=uigetdir % the folder to open is the folder with the soundsfiles. uigetdir chooses the directory

oldwav=[];

for N=1:length(dataX(:,1)) % This specifies that the code will run for the length of datax columns
    if N<=length(filelist)-1 % this is the part of the code making shure that everything under here will repeat itself for each row in the dataset
    
    end
    wavname= convertStringsToChars(filelist(N)); % convert strings (1) to characters(matlab counts each letter)
    loadwav=sprintf('%s.wav',wavname(10:end-4)); % these two lines linkes the wav.files with the file column. as the names matches we just take away what doesnt match. 
    TF = strcmp(loadwav,oldwav); % the TF just keeps the loop from loading if it recognises the file as the same så it doesnt load all over again everytime it changes row.
    if TF==0 % if its the same file, do this next section for the next row
        [data,fs]=audioread(fullfile(selpath,loadwav)); %loads the connected wav file
        %data2 = highpass(data(:,channel(N)),500,fs); % puts a highpass on  it (just necssesary for clicks)
        samp2=round(dataX*fs); % uses the start/endtime to just get the sections I want
    end
    oldwav=loadwav;
    datasample=data(samp2(N,1)-10:samp2(N,2)+10); % this adds 10 samples on each side of the cut to make sure everything is within. 
    
end 

I know that the code works in finding the right sections and extracting the data i want. What I want to do is to include a code in the loop which extracts the relative amplitude within each frequency from each of the signals into a matrix which I will use to get the mean relative power (y) which i will use for the powerspectrumplot below ( I already have x (0-48000 Hz):

figure 
xkhz = x/1000;
x2 = [xkhz, fliplr(xkhz)];
inBetween = [curve1, fliplr(curve2)];
fill(x2, inBetween, [.5 .5 .5]); % fills in the SD.
hold on;
plot(xkhz, y, 'r', 'LineWidth', 2);
xlim([0 20])
ylim([0 0.00000015])
xlabel ' Frequency (kHz)'
ylabel 'Relative power'
title ' Power Spectrum 100 random clicks'

I hope this was understandable, ask if you need more information. Thank you :)

Answer 1

[Mathieu NOE]

hello Anika

this is a starter... I modified a bit your code so I could easily make it work on my side . Simply comment the lines where you see the comment : % me and uncomment your original lines (that I commented).

So basically I copied some wav files in a new folder (anika) and run this code ; you get now an array of fft spectrum in array spectrum_all

like this (NB : amplitude in dB and freq in log scale)

now I am not understanding what you want to do next - the "relative power" code is not clear to me - can you expand on that please ?

my code (so far) :

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% me : generation of filelist from wav files list (d)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d = dir('*.wav');
N = numel(d);
for k = 1:N
    filelist(k,1) = string(d(k).name);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000;    % for maximal resolution take NFFT = signal length (but then no averaging possible)
OVERLAP = 0.75; % overlap will be used only if NFFT is shorter than signal length


% selpath=uigetdir % the folder to open is the folder with the soundsfiles. uigetdir chooses the directory
selpath=pwd; % me

oldwav=[];
for k=1:N % me %  you => % length(dataX(:,1)) % This specifies that the code will run for the length of datax columns
    
%     if k<=length(filelist)-1 % this is the part of the code making shure that everything under here will repeat itself for each row in the dataset
%     
%     end
    
    wavname= convertStringsToChars(filelist(k)); % convert strings (1) to characters(matlab counts each letter)
    
%     loadwav=sprintf('%s.wav',wavname(10:end-4)); % these two lines linkes the wav.files with the file column. as the names matches we just take away what doesnt match. 
    loadwav=sprintf('%s.wav',wavname(1:end-4)) % me 

    TF = strcmp(loadwav,oldwav); % the TF just keeps the loop from loading if it recognises the file as the same så it doesnt load all over again everytime it changes row.
    if TF==0 % if its the same file, do this next section for the next row
        [data,fs]=audioread(fullfile(selpath,loadwav)); %loads the connected wav file
        %data2 = highpass(data(:,channel(N)),500,fs); % puts a highpass on  it (just necssesary for clicks)
        
%        samp2=round(dataX*fs); % uses the start/endtime to just get the sections I want
        samp2(k,1) = round(numel(data)*0.3); % me
        samp2(k,2) = round(numel(data)*0.7); % me
        
    end
    oldwav=loadwav;
    datasample=data(samp2(k,1)-10:samp2(k,2)+10); % this adds 10 samples on each side of the cut to make sure everything is within. 
    
    %% now the averaged FFT spectrum code
    [freq, spectrum_raw] = myfft_peak(datasample,fs,NFFT,OVERLAP);
    ind = freq>0;
    freq = freq(ind);
    spectrum_all(:,k) = spectrum_raw(ind); % create an array of spectra (column oriented)
    leg{k} = loadwav; % create one legend per wav file loaded
    
end 

%% plot
    figure,semilogx(freq,20*log10(spectrum_all));grid on
    df = freq(2)-freq(1); % frequency resolution 
    title(['Averaged FFT Spectrum  / Fs = ' num2str(fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
    xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
    legend(leg);

%% functions are below this line (do not modify)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    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*ones(1,channels));
        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

[Mathieu NOE]

Hello @Annika

problem solved ?