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

16 views (last 30 days)
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 :)

Answers (1)

Mathieu NOE
Mathieu NOE on 25 Jan 2024
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

Community Treasure Hunt

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

Start Hunting!