Extract the relative amplitude from various signals to a matrix for a powerspectrum plot
15 views (last 30 days)
Show older comments
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 :)
0 Comments
Answers (1)
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
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!