Test-case: Add to check-alsabat.sh additional audio test

This patch adds run of check_wav_file_snr.m script for alsabat
capture file to analyze with STFT spectra the audio quality.
Dropped signal-to-noise ratio in a single FFT indicates presence
of glitch. A low SNR result triggers run of additional glitch
periodicity finder.

Signed-off-by: Seppo Ingalsuo <[email protected]>

test-case/check-alsabat.sh

@@ -19,7 +19,8 @@
 rm -f /tmp/bat.wav.*
 
 # shellcheck source=case-lib/lib.sh
-source "$(dirname "${BASH_SOURCE[0]}")"/../case-lib/lib.sh
+TESTDIR=$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)
+source "$TESTDIR/case-lib/lib.sh"
 
 OPT_NAME['p']='pcm_p'     	OPT_DESC['p']='pcm for playback. Example: hw:0,0'
 OPT_HAS_ARG['p']=1          	OPT_VAL['p']=''
@@ -95,6 +96,22 @@ function __upload_wav_file
     done
 }
 
+function check_wav_file_snr
+{
+    for file in /tmp/bat.wav.*
+    do
+	if test -s "$file"; then
+	    dlogi "Checking wav file $file"
+	    cd "$TESTDIR"/tools || {
+		die "Failed to enter directory tools"
+	    }
+	    octave --silent --no-gui --eval "check_wav_file_snr('$file', $frequency, '$LOG_ROOT/')" || {
+		die "Error: Script check_wav_file_snr.m found issues in $file"
+	    }
+	fi
+    done
+}
+
 # check the PCMs before alsabat test
 dlogi "check the PCMs before alsabat test"
 aplay   -Dplug$pcm_p -d 1 /dev/zero -q || die "Failed to play on PCM: $pcm_p"
@@ -119,4 +136,6 @@ alsabat -C$pcm_c -c $channel_c -r $rate -f $format -F $frequency -k $sigmak || {
         exit 1
 }
 
+check_wav_file_snr
+
 wait $playPID

tools/check_wav_file_snr.m

@@ -0,0 +1,355 @@
+function pass = check_wav_file_snr(fn, test_tone_f, logdir)
+
+if nargin < 3
+	logdir = '/tmp';
+end
+
+if nargin < 2
+	% Automatic find of tone frequency
+	test_tone_f = 0;
+end
+
+if nargin < 1
+	error('Need to provide file name to analyze');
+end
+
+if exist('OCTAVE_VERSION', 'builtin')
+	pkg load signal;
+end
+
+%fp = '/home/singalsu/Downloads/';
+%wf = 'bat.wav.3Vbz0k'; % total mess
+%wf = 'bat.wav.zhJhMf'; % glitch in right channel
+%wf = 'bat.wav.Wodg37'; % low quality
+%wf = 'bat.wav.WR5of5'; % Low signal level but OK
+%wf = 'bat.wav.05IgGp'; % good
+%wf = 'bat.wav.3UiPoZ'; % good
+%fn = fullfile(fp, wf);
+
+param.t_ignore_from_start = 30e-3;
+param.t_ignore_from_end = 0;
+param.t_step = 1.0e-3;
+param.n_fft = 1024;
+param.max_snr_drop = 6;
+param.min_snr = 50;
+test_passed = true;
+do_plot = 1;
+do_print = 1;
+nfig = 0;
+
+% Read audio
+fprintf(1, 'Loading file %s\n', fn);
+[x, fs, frames, channels] = read_audio(fn);
+[~, basefn, extfn] = fileparts(fn);
+
+for ch = 1 : channels
+	% Sanity check
+	if ~signal_found(x(:, ch))
+		continue;
+	end
+
+	% STFT
+	[stft, f, t, param] = compute_stft(x(:, ch), fs, param);
+
+	% Plot
+	if do_plot
+		fnstr = sprintf('%s%s', basefn, extfn);
+		idstr = sprintf('%s ch%d', fnstr, ch);
+		nfig = nfig + 1;
+		plot_specgram(stft, f, t, idstr);
+	end
+
+	% Check frequency, get STFT frequency index
+	[signal_idx, tone_f] = find_test_tone(stft, param, f, t, test_tone_f);
+
+	% Get levels and SNR estimates
+	meas = get_tone_levels(stft, param, f, t, signal_idx);
+
+	% Plot SNR
+	if do_plot
+		plot_levels(meas, t, idstr);
+		nfig = nfig + 1;
+	end
+
+	% Checks for levels data
+	meas = check_tone_levels(stft, param, meas, t);
+
+	% If poor SNR check if there's periodic glitches
+	meas = check_glitch_periodicity(x(:,ch), param, meas, ch);
+
+	% Plot glitch waveform
+	if do_plot && meas.num_glitches > 0
+		plot_glitch(x(:, ch), param, meas, idstr);
+		nfig = nfig + 1;
+	end
+
+	all_meas(ch) = meas;
+	if meas.success == false
+		test_passed = false;
+	end
+end
+
+fprintf(1, 'Ch, Pass, Frequency, SNR min, SNR avg, Signal avg, Noise avg, Noise max, Num glitch, 1st glitch\n');
+
+for ch = 1 : channels
+	meas = all_meas(ch);
+	fprintf(1, '%d, %2d,     %5.0f,  %6.1f,  %6.1f,     %6.1f,    %6.1f,    %6.1f,     %6d,  %8.3f\n', ...
+		ch, meas.success, tone_f, meas.min_snr_db, meas.mean_snr_db, ...
+		meas.mean_signal_db, meas.mean_noise_db, ...
+		meas.max_noise_db, meas.num_glitches, meas.t_glitch);
+end
+
+if do_plot && do_print && ~test_passed
+	for i = 1 : nfig
+		pfn = sprintf('%s/analyze_%s_%d.png', logdir, fnstr, i);
+		figure(i);
+		print(pfn, '-dpng');
+	end
+end
+
+if test_passed
+	fprintf(1, 'Passed\n');
+else
+	error('Failed');
+end
+
+end
+
+%
+% Helper functions
+%
+
+function meas = check_glitch_periodicity(x, param, meas, ch)
+
+if meas.success
+	return
+end
+
+i1 = round(param.fs * param.t_ignore_from_start);
+i2 = length(x) - round(param.fs * param.t_ignore_from_end);
+[b, a] = butter(2, 0.95, 'high');
+y = abs(filter(b, a, x));
+z = y(i1:i2);
+thr = mean(z) + std(z);
+[~, locs]= findpeaks(z, 'MinPeakHeight', thr);
+dlocs = diff(locs);
+if isempty(dlocs)
+	return
+end
+
+[counts, centers] = hist(dlocs);
+[counts, i] = sort(counts, 'descend');
+centers = centers(i) * 1e3 / param.fs;
+thr = mean(counts) + std(counts);
+idx = find(counts > thr);
+if ~isempty(idx)
+	fprintf(1, 'Ch%d periodic glitches possibly every', ch)
+	for i = 1 : length(idx)
+		fprintf(' %.1f ms (n = %d)', centers(idx(i)), counts(idx(i)));
+	end
+	fprintf(1, '\n');
+end
+
+end
+
+function success = signal_found(x)
+
+% All zeros or DC
+if abs(min(x) - max(x)) < eps
+	success = 0;
+else
+	success = 1;
+end
+
+end
+
+function meas = get_tone_levels(stft, param, f, t, signal_idx)
+
+signal_i1 = signal_idx - param.win_spread;
+signal_i2 = signal_idx + param.win_spread;
+if signal_i1 < 1
+	error('Too low tone frequency, increase FFT length');
+end
+
+signal_db = zeros(param.n_stft, 1);
+noise_db = zeros(param.n_stft, 1);
+snr_db = zeros(param.n_stft, 1);
+for i = 1 : param.n_stft
+	% Integrate signal power
+	p_signal = sum(stft(signal_i1 : signal_i2, i));
+
+	% Integrate noise power, but replace DC and signal with
+	% average noise level.
+	noise = stft(:, i);
+	noise_avg = mean(noise(signal_i2 : end));
+	noise(1 : param.win_spread) = noise_avg;
+	noise(signal_i1 : signal_i2) = noise_avg;
+	p_noise = sum(noise);
+
+	% Sign, noise, and "SNR" as dB
+	signal_db(i) = 10*log10(p_signal);
+	noise_db(i) = 10*log10(p_noise);
+	snr_db(i) = signal_db(i) - noise_db(i);
+end
+
+meas.noise_db = noise_db - param.win_gain;
+meas.signal_db = signal_db - param.win_gain;
+meas.snr_db = signal_db - noise_db;
+
+i1 = find(t > param.t_ignore_from_start, 1, 'first');
+i2 = find(t <= t(end) - param.t_ignore_from_end, 1, 'last');
+
+meas.mean_signal_db = mean(meas.signal_db(i1 :i2));
+meas.mean_noise_db = mean(meas.noise_db(i1 :i2));
+meas.mean_snr_db = mean(meas.snr_db(i1 :i2));
+meas.max_noise_db = max(meas.noise_db(i1 :i2));
+meas.min_snr_db = min(meas.snr_db(i1 :i2));
+
+end
+
+function meas = check_tone_levels(stft, param, meas, t)
+
+meas.t_glitch = 0;
+meas.num_glitches = 0;
+meas.success = true;
+
+% Find glitches from SNR curve drops
+i1 = find(t > param.t_ignore_from_start, 1, 'first');
+i2 = find(t <= t(end) - param.t_ignore_from_end, 1, 'last');
+idx = find(meas.snr_db(i1:i2) < meas.mean_snr_db - param.max_snr_drop);
+
+if ~isempty(idx)
+	idx = idx + i1 - 1;
+	didx = diff(idx);
+	meas.num_glitches = 1 + length(find(didx > 2));
+	start_idx = idx(1);
+	cidx = find(didx(1:end) > 1, 1, 'first');
+	if isempty(cidx)
+		end_idx = idx(end);
+	else
+		end_idx = idx(cidx);
+	end
+	meas.t_glitch = param.t_step * mean([start_idx end_idx] - 1) + ...
+		0.5 * param.n_fft / param.fs;
+	meas.success = false;
+end
+
+if meas.min_snr_db < param.min_snr
+	meas.success = false;
+end
+
+end
+
+function [signal_idx, tone_f] = find_test_tone(stft, param, f, t, test_tone_f)
+
+if test_tone_f > 0
+	err_ms = (f - test_tone_f) .^ 2;
+	signal_idx = find(err_ms == min(err_ms));
+	tone_f = f(signal_idx);
+	return
+end
+
+i1 = find(t > param.t_ignore_from_start, 1, 'first');
+i2 = find(t <= t(end) - param.t_ignore_from_end, 1, 'last');
+
+signal_idx_all = zeros(i2 - i1 + 1, 1);
+for i = i1 : i2
+	signal_idx_all(i - i1 + 1) = find(stft(:, i) == max(stft(:, i)),1, 'first');
+end
+
+signal_idx = round(mean(signal_idx_all));
+tone_f = f(signal_idx);
+
+end
+
+function [x, fs, frames, channels] = read_audio(fn)
+
+[x, fs] = audioread(fn);
+sx = size(x);
+frames = sx(1);
+channels = sx(2);
+
+end
+
+function [stft, f, t, param] = compute_stft(x, fs, param)
+
+sx = size(x);
+if sx(2) > 1
+	error('One channel only');
+end
+
+frames = sx(1);
+win = kaiser(param.n_fft, 20);
+param.win_spread = 7;
+param.win_gain = -13.0379;
+param.fs = fs;
+
+param.n_step = fs * param.t_step;
+param.n_stft = floor((frames - param.n_fft) / param.n_step);
+n_half_fft = param.n_fft / 2 + 1;
+scale = 1 / param.n_fft;
+f = linspace(0, fs/2, n_half_fft);
+t = (0 : (param.n_stft - 1)) * param.t_step;
+stft = zeros(n_half_fft, param.n_stft);
+
+for i = 1 : param.n_stft
+	i1 = (i - 1) * param.n_step + 1;
+	i2 = i1 + param.n_fft - 1;
+	s1 = fft(x(i1 : i2) .* win) * scale;
+	s2 = s1(1 : n_half_fft);
+	s = s2 .* conj(s2);
+	stft(:, i) = s;
+end
+
+end
+
+function fh = plot_specgram(stft, f, t, tstr)
+
+fh = figure;
+clims = [-160 0];
+imagesc(1e3 * t, f, 10*log10(stft), clims)
+set(gca, 'ydir', 'normal');
+colorbar;
+grid on;
+title(tstr, 'interpreter', 'none');
+xlabel('Time (ms)');
+ylabel('Frequency (Hz)');
+
+end
+
+function fh = plot_levels(meas, t, idstr)
+
+t_ms = 1e3 * t;
+fh = figure;
+subplot(3, 1, 1);
+plot(t_ms, meas.snr_db);
+grid on;
+ylabel('SNR (dB)');
+title(idstr, 'interpreter', 'none');
+subplot(3, 1, 2);
+plot(t_ms, meas.signal_db);
+grid on;
+ylabel('Signal (dBFS)');
+subplot(3, 1, 3);
+plot(t_ms, meas.noise_db);
+grid on;
+ylabel('Noise (dBFS)');
+xlabel('Time (ms)');
+
+end
+
+function fh = plot_glitch(x, param, meas, idstr)
+
+fh = figure;
+t_ms = 1e3 * (0:(length(x) - 1)) / param.fs;
+plot(t_ms, x);
+t_start = meas.t_glitch - 2 * param.t_step;
+t_end = meas.t_glitch + 2 * param.t_step;
+ax = axis();
+axis([1e3 * t_start 1e3 * t_end ax(3:4)]);
+grid on;
+title(idstr, 'interpreter', 'none');
+xlabel('Time (ms)');
+ylabel('PCM sample values');
+
+end