calibration tools added

subtools/calibration/TDT-analysescripts/TDT_Calibrator_Quiet_Analysis_Script.m

+% Below script was used to:
+% - Analyse calibration pulses
+% - Calculate background noise level
+% - Generate plots that show the calibration pulse and background noise
+% level
+
+% Written by: JB van der Heijdt
+% Date: April 2019
+
+clearvars;
+% Calibration
+CAL = load('JH-2019-03-21-calibrator.mat');
+
+cal   = CAL.Sounds.recorded; % [V]
+% [B,A]	= butter(4,60/CAL.fsRZ6,'high'); % Filtering removes about 0.5% from
+% overall rms. Makes sense based on 40 dB difference of LF components
+% cal		= filtfilt(B,A,cal);
+
+Ppulse  = rms(cal);
+DBref = 94 - 0.3; % dB SPL
+Pref  = 10^((DBref)/20); % [Pa] -0.3 dB correction for free-field microphones
+
+scaling  = Pref/Ppulse;
+Pcal = scaling*cal;
+
+nfftCal = 2^(nextpow2(length(cal)-1));
+
+[Farray,Fcal]	=	genfft(Pcal,nfftCal,CAL.fs);
+Fcal			=	sqrt(2)/length(cal) * Fcal;
+Fdb				=	20*log10(Fcal/P0);
+
+%% PLOT: Calibrator pulse and scaling
+fs=CAL.fs;
+hFig = figure(3);
+ax = subplot(1,3,1);
+plot(1/fs:1/fs:12,cal);
+xlabel('Time (s)','FontSize',24);
+ylabel('Amplitude','FontSize',24);
+title('Calibrator pulse','FontSize',24);
+ax.FontSize = 24;
+axis square
+
+ax = subplot(1,3,2);
+plot(1/fs:1/fs:0.012,cal(1:length(cal)/1e3));
+xlabel('Time (s)','FontSize',24);
+ylabel('Amplitude','FontSize',24);
+title('Calibrator pulse (detail)','FontSize',24);
+axis square
+ax.FontSize = 24;
+ax = subplot(1,3,3);
+[f,a] = getpower(Pcal,fs,'display',false,'nfft',nfftCal);
+semilogx(f,a,'color',[0 0 0]);
+grid on
+xlabel('Frequency (kHz)','FontSize',24);
+ylabel('Sound level (dB SPL)','FontSize',24);
+title('Power spectrum','FontSize',24);
+axis square
+ax.FontSize = 20;
+		set(gca,'XTick',[0.05 1  5 10]*1000,...
+			'XTickLabel',[0.05 1  5 10]);
+w = aweight(f);
+m = a+w;
+		
+p		= 10.^(a/20);
+psum	= sqrt(sum(p.^2));
+L		= 20*log10(psum);
+
+p		= 10.^(m/20);
+psum	= sqrt(sum(p.^2));
+La		= 20*log10(psum);
+
+levels = [20*log10(rms(Pcal)) L La]
+
+str = ['L_{rms} = ' num2str((levels(1)))];
+text(1000+100,94-3,str,'FontSize',20);
+
+% nicegraph(gcf)
+
+%% Analysis: quiet (has different scaling value, too)
+CAL = load('JH-linsweep-2019-03-28-calibrator.mat');
+cal   = CAL.Sounds.recorded; % [V]
+
+Ppulse  = rms(cal);
+P0 = 1;
+DBref = 94 - 0.3; % dB SPL
+Pref  = P0*10^((DBref)/20); % [Pa] -0.3 dB correction for free-field microphones
+
+scaling  = Pref/Ppulse;
+Pcal = scaling*cal;
+
+% nfftCal = 2^(nextpow2(length(cal)-1));
+%  for ii = 1:19
+% 	[Farray,Fcal]	=	genfft(Pcal(floor((ii:ii+1)*length(Pcal)/20)),nfftCal,CAL.fsRZ6);
+% 	Fcal	=	sqrt(2)/length(cal) * Fcal;
+% 	Fdb(ii,:)			= 20*log10(Fcal/P0);
+%  end
+%
+ % Load background noise recording
+load('/Users/jasperhvanderheijdt/Documents/Werk/Meetdata/Sweeps/JH-9998-2019-03-27/JH-linsweep-8888-2019-03-27-0001-quiet.mat')
+
+recordQuiet = scaling*Sounds.recorded;
+recordQuietA = filterA(recordQuiet,fsRZ6);
+% [b,a]=butter(4,100/fsRZ6,'high');
+% recordQuietFilt=filter(b,a,recordQuiet);
+overallNoise = 20*log10(rms(recordQuietA))
+%
+nfft = 2^(nextpow2(length(recordQuiet))-1);
+
+[fArray,fQuiet] = genfft(recordQuiet,nfft,fsRZ6);
+Nlength = length(fArray)/2 + 1; % remove mirrored FFT
+fArray = fArray(1:Nlength);
+fQuiet = fQuiet(1:Nlength);
+fQuiet = sqrt(2)/length(recordQuiet) * fQuiet; % Correction (Parseval theorem)
+fQuiet= 20*log10(fQuiet/P0); % dB conversion
+tArray = (1:length(recordQuiet))/fsRZ6;
+
+% Now use Welch' method (windowing reduces noise)
+nWin = 300; % number of windows (300 windows = 0.1 s for 30s recording).
+win = length(recordQuiet)/nWin;
+nfft=[];
+[pRaw,freqArray] = pwelch(recordQuiet,win,[],nfft,fsRZ6);
+pQuiet = 20*log10(pRaw);
+% [pFilt,~] = pwelch(recordQuietFilt,win,[],nfft,fsRZ6);
+% pFilt = 20*log10(pFilt);
+
+
+% Add a-weighting filter
+[W,A]=aweight(freqArray); % A-weighted dB value (additive)
+pQuietA = W'+pQuiet;
+
+% preset labels etc.
+xVal = [1e2 5e2 1e3 5e3 1e4];
+xName = [0.1 0.5 1 5 10]; % kHz
+xLabName = 'Frequency (kHz)';
+
+% Plotting
+hFig2 = figure(2);
+clf;
+hAx(1) = subplot(1,3,1); % Time domain unscaled
+plot(tArray,recordQuiet/scaling);
+xlabel('Time (s)','FontSize',24);
+ylabel('Unscaled mplitude (V)','FontSize',24);
+title('Time domain','FontSize',24);
+
+hAx(2) = subplot(1,3,2); % dB SPL
+semilogx(freqArray,pQuiet); 
+% semilogx(fArray,fQuiet);
+xlim([10 fsRZ6/2]);
+xlabel(xLabName,'FontSize',24);
+ylabel('Sound level (dB SPL)','FontSize',24);
+title('Frequency Domain','FontSize',24);
+xticks(xVal);
+xticklabels(xName);
+
+hAx(3) = subplot(1,3,3); % A-weighted
+semilogx(freqArray,pQuietA); 
+% semilogx(freqArray,pFilt);
+hold on
+% semilogx(fArray,fQuiet+W);
+
+xlim([10 fsRZ6/2]);
+xlabel(xLabName,'FontSize',24);
+ylabel('Sound level (dB A)','FontSize',24);
+title('A-weighted frequency domain','FontSize',24);
+xticks(xVal);
+xticklabels(xName);
+
+% Xticks fontsize
+for iAx = 1:length(hAx)
+	hAx(iAx).FontSize = 24;
+end
+hTxt = suptitle('Recording of background noise');
+hTxt.FontSize = 36;
+str = ['L_{rms} = ' num2str(round(overallNoise,1)) ' dB A'];
+text(15,0,str,'FontSize',24);
+fancify(hFig2)
+

subtools/calibration/TDT-analysescripts/TDT_GWN_Analysis_Script.m

+CAL = load('/Users/jasperhvanderheijdt/Documents/Werk/Meetdata/Sweeps/JH-9998-2019-03-27/JH-linsweep-8888-2019-03-27-0001-calibrator.mat');
+CAL = load('load('JH-linsweep-2019-03-28-calibrator.mat');
+cal   = CAL.Sounds.recorded; % [V]
+
+Ppulse  = rms(cal);
+P0 = 1;
+DBref = 94 - 0.3; % dB SPL
+Pref  = P0*10^((DBref)/20); % [Pa] -0.3 dB correction for free-field microphones
+
+scaling  = Pref/Ppulse;
+Pcal = scaling*cal;
+
+nfftCal = 2^(nextpow2(length(cal)-1));
+
+[Farray,Fcal]	=	genfft(Pcal,nfftCal,CAL.fsRZ6);
+Fcal	=	sqrt(2)/length(cal) * Fcal;
+Fdb			= 20*log10(Fcal/P0);
+
+%% Some sound recording 
+load('JH-0000-19-03-27-0000-0003.sphere','-mat');
+
+% Scaling, FFT's, corrections
+
+nfft = 2^(nextpow2(length(soundRecording))-1);
+soundRecording = scaling*soundRecording;
+
+[fArray,fBurst] = genfft(soundRecording,nfft,cfg.RZ6Fs);
+fBurst = sqrt(2)/length(soundRecording) * fBurst; % Correction (Parseval theorem)
+fBurst= 20*log10(fBurst/P0); % dB conversion
+
+figure(1);
+clf;
+semilogx(fArray,fBurst);
\ No newline at end of file

subtools/calibration/TDT-analysescripts/TDT_ImpulseResponse_Analysis_Script.m

+% Analysis of Impulse responses of spherelab
+% TODO: 
+% - Measure theoretical distance of initial response and first
+% - Calibrate soundlevels (does this work in time-domain?)
+% reflection (with ruler, verify technique)
+% - Determine size and timing of initial response of 1st reflection in
+% measurements. Manually or automated (peak detection algorithm?)
+% Pictures: pulse provided. Typical response
+
+
+% First things first: scaling factor (copye-pastable)
+fdir = '/Users/jasperhvanderheijdt/Documents/Werk/Meetdata/Click/2019-04-05';
+
+
+CAL = load('JH-2019-04-05-calibrator.mat');
+% CAL   = load('/Users/jasperhvanderheijdt/Documents/Werk/Meetdata/Aco-pacific test/ACOpacific_calibration.mat');
+cal   = CAL.Sounds.recorded; % [V]
+% cal		= CAL.RZ6_mic(:,1);
+
+% [B,A]	= butter(4,60/CAL.fsRZ6,'high'); % Filtering removes about 0.5% from
+% overall rms. Makes sense based on 40 dB difference of LF components
+% cal		= filtfilt(B,A,cal);
+
+Ppulse  = rms(cal);
+% P0 = 2e-5; % [dB] 20 uPa!
+ P0 = 1;
+DBref = 94 - 0.3; % dB SPL
+Pref  = P0*10^((DBref)/20); % [Pa] -0.3 dB correction for free-field microphones
+
+scaling  = Pref/Ppulse;
+Pcal = scaling*cal;
+
+nfftCal = 2^(nextpow2(length(cal)-1));
+
+[Farray,Fcal]	=	genfft(Pcal,nfftCal,CAL.fsRZ6);
+Fcal	=	sqrt(2)/length(cal) * Fcal;
+Fdb		= 20*log10(Fcal/P0);
+
+%% Get speakerdistances
+spkDistance = importfile('speakerdistance.csv',2,30);
+% Calculate expected delays
+v=343; % [m/s]
+
+for ii = 1:29
+	x= spkDistance.distspkmic(ii); % [m]
+	y = spkDistance.distspwall(ii); % [m]
+	delayIR(ii) = 1e3*x/v; % [ms]
+	delayReflection(ii) = 1e3*(x+2*y)/v; % [ms]
+end
+
+%% Plot example time domain
+ % with inputs
+fdir = '/Users/jasperhvanderheijdt/Documents/Werk/Meetdata/Click/2019-04-05/';
+fnameBase = 'JH-2019-04-05-click-Spk';
+fExt	= '.mat';
+
+% test: load 1 file
+iSp = 6;
+filename = [fdir fnameBase num2str(iSp) fExt];
+load(filename);
+
+% construct time Array
+tArray = 1e3*(1:length(Sounds.recorded))/fs; % [ms]
+tWin = (tArray > 100 & tArray < 120); % window of interest for peak detection
+
+
+% remove all peaks less than 'x'% of max value
+% highest = max(peaks);
+% xFactor = .1;
+% dum = peaks > xFactor*highest;
+% peaks = peaks(dum);
+% locs = locs(dum);
+
+% some plotting parameters
+xText = 'Time (ms)';
+yText = 'Response';
+xlimits = [90 110];
+ylimits = [-2 2];
+txtBase0 = 'Pulse start: ';
+txtBase1 = 'recorded delay IR: ';
+txtBase2 = 'recorded delay 1st reflection: ';
+txtBase3 = 'theoretical delay IR: ';
+txtBase4 = 'theoretical delay 1st reflection: ';
+xstartTxt = 90;
+prec = 2; % decimals - rounding precision
+hFig = figure(1);
+for iSp = 1:29
+	filename = [fdir fnameBase num2str(iSp) fExt];
+	load(filename);
+	
+	% peak detection
+	rec = Sounds.recorded;
+	[peaks,locs] = findpeaks(rec(tWin),tArray(tWin),'MinPeakProminence',0.4,'MinPeakHeight',1e-3);
+	[peakPulse,locPulse] = findpeaks(Sounds.out(tWin),tArray(tWin),'MinPeakProminence',15,'MinPeakHeight',1e-3);
+	subplot(3,1,1);
+	plot(tArray,Sounds.in);
+	xlim(xlimits);
+	xlabel(xText);
+	ylabel(yText);
+	title('Input: 1 ms pulse');
+	subplot(3,1,2);
+	hold on
+	plot(tArray,Sounds.out);
+	plot(locPulse(1),peakPulse(1),'o');
+	xlim(xlimits);
+	title('Signal sent to speaker');
+	xlabel(xText);
+	ylabel(yText);
+	text(xstartTxt,8,[txtBase0 num2str(round(locPulse(1),1)) ' ms'],'FontSize', 20);
+	hold off
+	subplot(3,1,3);
+	plot(tArray,rec);
+	hold on;
+	plot(locs,peaks,'o');
+	title('Recorded response');
+	xlabel(xText);
+	ylabel(yText);
+	xlim(xlimits);
+	ylim(ylimits);
+	text(xstartTxt,1,   [txtBase1 num2str(round(locs(1)-locPulse(1),prec)) ' ms'],'FontSize', 20);
+	text(xstartTxt,0.5, [txtBase2 num2str(round(locs(2)-locPulse(1),prec)) ' ms'],'FontSize',20);
+	text(xstartTxt,-0.5,[txtBase3 num2str(round(delayIR(iSp),prec)) ' ms'],'FontSize', 20);
+	text(xstartTxt,-1,  [txtBase4 num2str(round(delayReflection(iSp),prec)) ' ms'],'FontSize',20);
+	hold off
+
+	fancify(hFig)
+	hTxt = suptitle(['Impulse Response of speaker #' num2str(iSp)]);
+	hTxt.FontSize = 36;
+
+	pause
+end

subtools/calibration/TDT-analysescripts/TDT_Noise_analysis_I.m

+% Script for analyzing the data obtained on 2019-April- 17th
+% Goal: troubleshooting.
+% Note: this is an ugly file with lots of dead ends and dummy parameters. Fix before presenting
+
+% Step 1: Calibrator pulse
+% Step 2: Load data files and separate per speaker and per frequency band
+% Step 3: Calculate scaled RMS values in dB SPL for each speaker and
+% frequency band
+% Plot the findings
+ % first set: example of all 3 Noise Bursts in T- and F-domain (any
+ % speaker), at each step of analysis
+ % 2nd set: X-axis: dB attenuation; Y-axis: RMS value (dB) for
+ % Sounds.recording - to check for any irregularities.
+ % 3rd set: if irregularities occur: plot Time- and F-domain 
+ 
+ %%%%%%%%%%%%%%%%%%%%%%%Let's get to work%%%%%%%%%%%%%%%%%%%
+ 
+% Calibrator pulse
+calFileName = 'JH-2019-04-17-calibrator.mat';
+CAL = load(calFileName);
+cal = CAL.Sounds.recorded;
+Ppulse  = rms(cal);
+
+DBref = 94 - 0.3; % dB SPL
+Pref  = 10^((DBref)/20); % [Pa] -0.3 dB correction for free-field microphones
+
+scaling  = Pref/Ppulse;
+Pcal = scaling*cal;
+
+nfftCal = 2^(nextpow2(length(cal)-1));
+
+[Farray,Fcal]	=	genfft(Pcal,nfftCal,CAL.fsRZ6);
+Fcal	=	sqrt(2)/length(cal) * Fcal;
+Fcal		= 20*log10(Fcal);
+
+%% Load files and separate as desired
+
+	fdir = '/Users/jasperhvanderheijdt/Documents/Werk/Meetdata/test Sndburst/JH-3630-19-04-17/trial/';
+	filelist = dir(fdir); 
+	dirFlags = [filelist.isdir]; 
+	filelist(dirFlags) = []; % remove directories (including hidden directories)
+% 	Select desired recording parameters
+	atten	= 0:1:80; % dB
+	spk		= 1:3;
+	
+ a=0; b=0; c=0; d=0; e=0; f=0; g=0; h=0; k=0;
+ dump = a;
+% 	for iTrial = 1
+	for iTrial = 1:length(filelist)
+		TRIAL = load(filelist(iTrial).name,'-mat');
+		
+		% Extract parameters
+			% All sounds. Speaker info. Parameters info (noise type)
+		noiseType			= 	TRIAL.trialsingle.stim.parameters;
+		speakerPos			=	TRIAL.trialsingle.stim.X; % azimuth position
+		% Data into 9 subgroups 
+		switch speakerPos
+			case 10
+				if noiseType == 1
+					a=a+1;
+					sndFilt1(a,:) = TRIAL.Sounds.sndFilt1;
+					sndFilt2(a,:) = TRIAL.Sounds.sndFilt2;
+					soundIn(a,:)  = TRIAL.Sounds.played;
+					soundRec(a,:) = TRIAL.Sounds.recording;
+				elseif noiseType ==2
+					b=b+1;
+					sndFilt12(b,:) = TRIAL.Sounds.sndFilt1;
+					sndFilt22(b,:) = TRIAL.Sounds.sndFilt2;
+					soundIn2(b,:)  = TRIAL.Sounds.played;
+					soundRec2(b,:) = TRIAL.Sounds.recording;
+				elseif noiseType == 3
+					c=c+1;
+					sndFilt13(c,:) = TRIAL.Sounds.sndFilt1;
+					sndFilt23(c,:) = TRIAL.Sounds.sndFilt2;
+					soundIn3(c,:)  = TRIAL.Sounds.played;
+					soundRec3(c,:) = TRIAL.Sounds.recording;
+% 				else
+% 					error('Invalid speakerPos/noiseType combination');
+				end
+			otherwise
+% 				dump = dump+1;
+% 			case -10
+% 			case 10
+		end
+% 		clear('TRIAL'); % cleanup
+	end
+		soundRec = soundRec*scaling;
+		soundRec2 = soundRec2*scaling;
+		soundRec3 = soundRec3*scaling;
+	%% Now let's plot RMS first 
+	
+	% Strong highpass filter
+	ord = 4; % filter order
+	W = 100/48848;
+	[b,a] = butter (ord, W, 'high'); 
+	soundDum = filter(b,a,soundRec);
+% 	soundDum=soundRec;
+	% Welch' method (windowing reduces noise)
+	nWin = 50; % number of windows (50 windows = 0.1 s for 5s recording).
+	win = length(soundDum)/nWin;
+	nfft=[];
+	for iTrial = 1:length(atten)
+% 		[pRec(iTrial,:),freqArray] = pwelch(soundDum(iTrial,:),win,[],nfft,CAL.fsRZ6);
+		[fArray,pRec(iTrial,:)] = genfft(soundDum(iTrial,:),512,CAL.fsRZ6/2);
+	end
+	% pRaw   = pRaw(15:end);
+
+	SlevelF = 20*log10(rms(pRec(:,1.:end)'));
+	SlevelT = 20*log10(rms(soundDum'));
+% 	pRecdB = 20*log10(pRec);
+figure(1);
+clf
+	plot(atten,SlevelT,'*'); hold on
+	plot(atten,SlevelF,'r*');
+	xlabel('dB attenuation','FontSize',24);
+	ylabel('RMS level (dB SPL)','FontSize',24);
+	title('BB noise RMS levels','FontSize',24);
+	legend('T-domain RMS','F-domain RMS');
+%% troubleshoot visualizations
+pRecdB = 20*log10(pRec);
+figure(1);
+for ii = 1:81
+	clf;
+	xlim([-80 80]);
+	semilogx(freqArray(10:end),pRecdB(ii,10:end));
+	hold on
+	line([250 250], ylim);
+	line([750 750],ylim);
+	line([2e3 2e3], ylim);
+	line([6e3 6e3], ylim);
+	title(atten(ii));
+	text(100,20,num2str(Slevel(ii))); 
+	pause;
+end
+	
+

subtools/calibration/TDT-analysescripts/TDT_Sines_Analysis_Script.m

+% This script:
+% Analyses the recording of various Sines generated by the speakers connected to RZ6 system of
+% the patient lab
+% This script:
+% -loads each recording
+% - groups them by speaker and attenuation level
+% - calculates the fft's (using genfft helper function)
+% - visualizes Gain the F-domain (using suptitle and fancify helper functions) and finally
+% - save the resulting figure as a .png in the designated directory.
+% Written by: JB van der Heijdt
+
+clearvars;
+
+% Different measurement dates (check that all fit!)
+% 2019-03-21: complete measurement, without extra amplifiers [obsolete]
+% 2019-04-25: complete measurement, with extra amplifiers [current]
+% NOTE: for this script, check the dates for the Master file, fdir and the
+% calibrator file
+
+% Load master file & get important parameters
+
+	fnameM = 'JH-2019-04-25-master.mat';
+	Master = load(fnameM);
+	% fBase = ... etc.
+	fBase = Master.fBase;
+	 
+ % Prepare directory list
+	fdir = '/Users/jasperhvanderheijdt/Documents/Werk/Meetdata/sineTests/2019-04-25';
+	filelist = dir(fdir); 
+	dirFlags = [filelist.isdir]; 
+	filelist(dirFlags) = []; % remove directories (including hidden directories)
+% 	Select desired recording parameters
+% 	spNo = 1; % speaker number
+	atten = 0:10:80; % dB
+
+% Calibration
+
+% CAL = load('sineTestsJH-2019-03-12-calibrator.mat');
+CAL = load('JH-2019-04-25-calibrator.mat');
+% CAL   = load('/Users/jasperhvanderheijdt/Documents/Werk/Meetdata/Aco-pacific test/ACOpacific_calibration.mat');
+cal   = CAL.Sounds.recorded; % [V]
+% cal		= CAL.RZ6_mic(:,1);
+
+% [B,A]	= butter(4,60/CAL.fsRZ6,'high'); % Filtering removes about 0.5% from
+% overall rms. Makes sense based on 40 dB difference of LF components
+% cal		= filtfilt(B,A,cal);
+
+Ppulse  = rms(cal);
+% P0 = 2e-5; % [dB] 20 uPa!
+ P0 = 1;
+DBref = 94 - 0.3; % dB SPL
+Pref  = P0*10^((DBref)/20); % [Pa] -0.3 dB correction for free-field microphones
+
+scaling  = Pref/Ppulse;
+Pcal = scaling*cal;
+
+nfftCal = 2^(nextpow2(length(cal)-1));
+
+[Farray,Fcal]	=	genfft(Pcal,nfftCal,CAL.fsRZ6);
+Fcal	=	sqrt(2)/length(cal) * Fcal;
+Fdb		= 20*log10(Fcal/P0);
+
+% Alternative: % [Fdb,a] = getpower(Pcal,fs,'display',true,'nfft',nfftCal);
+
+%% Data Extraction, FFT calculation, Visualization
+% For each speaker and each attenuation level (29x9=261) graphs
+% all frequencies at that level (7 subplots) in a 2x4 figure.
+
+% Set subplot parameters 
+	xLabName = 'Frequency (kHz)';
+	yLabName = 'Magnitude (dB SPL)';
+	titleExt	 = ' dB atten';
+	xlimits = [.1 1.9e4]; % [Hz]
+	GooderXTickLabel = [1 10 20];
+	xVals = [1000 10000 20000];
+% Set figure-specific parameters (name, supertitle)
+
+
+% Run the figure: Calculate FFT and plot 
+	load(filelist(1).name);
+	Sines = nan(7,length(Sounds(1).recorded)); % pre-allocate
+	nfft = 2^nextpow2(length(Sines(1,:)));
+
+% for spNo = 1:29
+for spNo = 1
+	hFig = figure(1);
+	%	Position parameters for suptitle (down below)
+	azPos = num2str(round(Master.speakers(spNo,4)));