Untitled

% Open wav files and read data for post-processing and plotting
% Use for Dyson and Donlim measured wav files from 10 mic chambers

%% [1] Select folder & subfolders containing .wav files
clearvars -except out
close all

cd(uigetdir);
wav_files = dir('**/*.wav');
wav_files=struct2table(wav_files);
wav_files=sortrows(wav_files,'datenum');

%% [2] File loop
for ii=1:height(wav_files)
    disp(['Computing file ',num2str(ii),' of ',num2str(height(wav_files))]);
    spectra=[];

    if height(wav_files) > 1
        filename=wav_files.name{ii};
        fldrname=wav_files.folder{ii};
    else
        filename=wav_files.name;
        fldrname=wav_files.folder;
    end
    dummy_txt=strsplit(filename,'\');
    measName=dummy_txt{end}; % assign legend entries
    % measName=filename;

    [dataAllChn,fs] = audioread(fullfile(fldrname,filename));

    % % remove empty channels
    dataAllChn(:,max(dataAllChn,[],1)<0.02)=[];

    if contains(filename,'calib')
        scalefactor = readNovaCalib(filename);
        scalefactor = repelem(scalefactor,1,size(dataAllChn,2));
    else
        % scalefactor = repelem(1,1,size(dataAllChn,2)); % use for calibration measurements

        % Donlim mic calibration
        cal = repelem(1,10); % no scaling
        % cal = repelem(72.2649,10); % 1000 Hz measured at Donlim X986 chamber with their calibrator
        cal = [74.6756,73.0778,72.5379,73.4923,74.828,72.8995,74.3011,73.5128,72.2771,71.9492]; % 1000 Hz measured at Donlim X938 chamber with B&K
        scalefactor = 10.^((94-cal)./20);
        scalefactor = scalefactor*1.5/2; % distance scaling from 1.5m radius hemisphere to match 2m data
    end

    dataAllChn=dataAllChn(:,1:10); % remove excess channels
    pScaled=dataAllChn.*scalefactor;

    if ~exist("out","var") % append new data from different folders
        out = struct;
        out(1).measName = measName;
        out(1).data = timetable(pScaled,SampleRate=fs);
    else
        out(end+1).measName = measName;
        out(end).data = timetable(pScaled,SampleRate=fs);
    end
end

%% plot data
t_roi = [5 10]; % time region of interest
df = 4;         % freq resolution for FFT and spectrogram
hemi_r = 2;     % hemisphere radius
measNameList={out.measName}.';

% select data points to plot
n = 1:numel(measNameList); % plot all
% n = find(contains(measNameList,'NoJug','IgnoreCase',true))'; % plot data points containing search token
% n=[1,2,12,13]; % plot specific data points

% harmonic magnitude extraction
H_order=1:7;    % harmonic orders to extract
Hx=zeros(numel(n),numel(H_order));
RPM = repelem(15000,numel(n),1); % use this if all wav files at 15 kRPM
% RPM=[25,50,75,100,125,150,175,200,225,250]*60; 

% Spectrogram and CPB parameters
overlap = 0.95;
fband=[20,25,31.5,40,50,63,80,100,125,160,200,250,315,400,500,650,800,1000,1250,1600,2000,2500,3150,4000,5000,6400,8000,10000,12500,16000,20000];

% % -----------------------------------------
close all
figure(1);clf
tiledlayout('flow','TileSpacing','tight')
figure(2);clf
tiledlayout('flow','TileSpacing','tight')
figure(3);clf

count=0;
for m=n
    count=count+1;
    fs = out(m).data.Properties.SampleRate;
    roi_idx = out(m).data.Time > seconds(t_roi(1)) & out(m).data.Time < seconds(t_roi(2));
    signal = out(m).data.pScaled(roi_idx,:);
    AwtFilt = weightingFilter('A-weighting',out(m).data.Properties.SampleRate);
    signalAwt = AwtFilt(signal);

    % % compute spectra
    [spectra,f]=frequencyspectrum(signal,fs,df,@hann,[],[],[],[]); % see MATE functions
    spectra=20*log10(spectra/20e-6); % convert to dB SPL
    spectra(:,end+1)=10.*log10(mean(10.^(spectra/10),2)); % append 10 channel average at end

    % % compute SWL from time series
    splNChn = 20*log10(sqrt(mean(signalAwt.^2))/2e-5); % compute SPL from time series
    splAvg = 10*log10(mean(10.^(splNChn/10)));  % hemispherical average SPL
    SWL(m) = round(splAvg + 10*log10(2*pi*hemi_r^2),1); % compute SWL from SPL

    % % channel loop
    for nn=1:size(signal,2)
        % compute channel CPB
        [cpbNChnT, fband, cpbT] = gene_ThirdOctave_levels(signalAwt(:,nn), fs, 0.25); % from MATE
        cpbNChn(:,nn) = mean(cpbNChnT,2); % use time-averaged CPB level

        % compute psychoacoustic metrics (takes long time!)
        % [psyML.Roughness(count,nn),psyML.Tonality(count,nn),psyML.Loudness(count,nn),psyML.Sharpness(count,nn),psyML.Fluc(count,nn),psyML.Annoy(count,nn),~] = PsychAc_MatlabMetrics(signal(:,nn),fs,1);
        % [psySQAT.Roughness(count,nn),psySQAT.Tonality(count,nn),psySQAT.Loudness(count,nn),psySQAT.Sharpness(count,nn),psySQAT.Fluc(m,nn),psySQAT.Annoy(count,nn),~] = PsychAc_SQATMetrics(signal(:,nn),fs,1);
    end
    
    % % compute Dyson psychoacoustic metric
    % psyDyson{count} = T11044_TonalLoudness_Nov2015('Pa',signal,'Fs',fs,'show',0);
    % Dyson_TLu(count)=max([psyDyson{count}.TL]);

    % % compute CPB and SWL
    CPB(:,m) = 10*log10(mean(10.^(cpbNChn/10),2)); % hemispherical average CPB
    % SWL(m) = 10*log10(sum(10.^(CPB(:,ii)/10))) + 10*log10(2*pi*hemi_r^2); % compute SWL from summing 1/3 octave bands   
    
    % % plot 10 mic spectra
    figure(1) 
    nexttile
    semilogx(f,spectra(:,1:end-1)); grid on; hold on
    xlim([100 20000])
    % ylim([-20 65])
    title(measNameList{m},'interpreter','none','FontSize',8)
    ylabel('SPL (dB re 20 \muPa)');xlabel('Frequency (Hz)')
    
    % % plot spectrogram for channel with highest SPL
    figure(2) 
    nexttile        
    window_size = 2^nextpow2(round((fs/df)*overlap/2)*2);
    windft = fft(hanning(window_size));
    fac = 2/sqrt(windft(1)^2);
    [~,NmaxSPL] = max(splNChn);
    [s, f, t, ps] = spectrogram(signal(:,NmaxSPL),hanning(window_size),floor(window_size*overlap),f,fs,'psd','yaxis');
    img2plot = 20*log10(abs(s)*fac/20e-6);
    imagesc(t,f(2:end),img2plot(2:end,:));
    set(gca,'YDir','Normal','YScale','linear','ylim',[4 10e3],'clim',[20 60])
    title(measNameList{m},'interpreter','none','FontSize',8)
    xlabel('Time (s)');ylabel('Frequency (Hz)')
    
    % % plot average spectrum
    figure(3)
    semilogx(f,spectra(:,end)); grid on; hold on
    xlim([100 10000])

    % % extract harmonic magnitudes
    Hf_win = [-10; 10] + (H_order * RPM(count)/60);
    Hf_win=floor(Hf_win/df);
    Hf_win(Hf_win<=0)=1;
    for mH = 1:numel(H_order)
        Hx(m,mH)=max(spectra(Hf_win(1,mH):Hf_win(2,mH),end));
    end
end

% create legend text
legendtxt = join([measNameList(n),string(SWL(n)'),repelem('dBA',numel(n),1)],' ');

% create truncated legend text 
if length(measNameList{end})>20
    legendtxtTrunc = join([extractBefore(measNameList(n),10),...
    repelem('...',numel(n),1),...
    extractAfter(measNameList(n),strlength(measNameList(n))-10),...
    string(SWL(n)'),repelem('dBA',numel(n),1)],' ');
else
    legendtxtTrunc = legendtxt;
end

figure(1); title('10 Mic Spectra')
figure(2); title('Spectrogram (Highest SPL Channel)')
figure(3); title('Average Spectrum')
legend(legendtxt,'interpreter','none')
ylabel('SPL (dB re 20 \muPa)');xlabel('Frequency (Hz)')

% % plot 1/3 octave bands
figure(4)
bar(CPB(:,n)); grid on
set(gca,'XTick',1:29,'XTickLabel',fband)
legend(legendtxt,'interpreter','none')
ylabel('SPL (dB re 20 \muPa)');xlabel('Frequency (Hz)')
title('1/3 Octave Bands')

% % plot harmonic magnitudes
figure(5)
Hx(sum(Hx,2)==0,:)=[];
plot(Hx,'Marker','o')
set(gca,'XTick',1:numel(n),'XTickLabel',legendtxtTrunc,'TickLabelInterpreter','none')
legend(strcat('H',string(H_order)'))
grid on
ylabel('Harmonic SPL (dBA)')
title('Harmonic Magnitudes')

% % plot TLu
% figure(6)
% bar(Dyson_TLu)
% set(gca,'XTick',1:numel(n),'XTickLabel',legendtxtTrunc,'TickLabelInterpreter','none')
% grid on
% ylabel('Dyson Tonal Loudness (TLu)')

%% Process psychoacoustic data
count=0;
for m=1:size(psychac_ML,1)
    % channel loop
    for nn=1:size(psychac_ML,2)
        % % % % % dummy(nn)=mean(psychac_ML(1:2,nn).loudness);
    end
end

%% edit spectrum and estimate SWL
spect_raw = out(12).data.spl_spect(end,:);
spect_smooth = spect_raw;

f_target = [500 750 1000 1250];
for n = 1:numel(f_target)
    win_ct = find(~floor((f-f_target(n))/10),1);
    win_width = 7;
    spect_win = smooth(spect_raw(win_ct-win_width:win_ct+win_width),...
        win_width*2,'sgolay',3);
    spect_smooth(win_ct-win_width:win_ct+win_width)=spect_win;
end

spectsum1 = 10*log10(sum(10.^(spect_raw/10))); % compute SPL from spectrum
spectsum2 = 10*log10(sum(10.^(spect_smooth/10))); % compute SPL from spectrum
del_spl = spectsum1 - spectsum2;

figure
semilogx(f,[spect_raw; spect_smooth])
xlim([min(f_target)-500 max(f_target)+500])
xlim([100 10000])
title("SWL reduction "+del_spl+" dB")
grid on

%% Extract order waveforms to analyze to Signal Analyzer App
m=2;
H=1:7;

t_roi=[1 2];
fs = out(m).data.Properties.SampleRate;
roi_idx = out(m).data.Time > seconds(t_roi(1)) & out(m).data.Time < seconds(t_roi(2));
signal = out(m).data(roi_idx,:);
signal = downsample(signal,10);
fs = fs/10;

xrec = orderwaveform(signal.pScaled(:,1),fs,repelem(15000,numel(signal(:,1))),H);

signalH = timetable(signal.Time,xrec);