%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sound Analysis using MATLAB %
% %
% Author: M.Sc. Eng. Hristo Zhivomirov 10/29/12 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear, clc, clf
% get a section of the sound file
[x,fs] = wavread('example.wav', [20000 65001]);
soundsc(x,fs);
x = x(:,1); % get the first channel
xmax = max(abs(x)); % find the maximum value
x = x/xmax; % scalling the signal
% time & discretisation parameters
N = length(x);
t=(0:N-1)/fs;
% plotting of the waveform
figure(1)
plot(t,x)
xlim([0 max(t)])
ylim([-1.1*max(abs(x)) 1.1*max(abs(x))])
grid on
set(gca, 'FontName', 'Times New Roman', 'FontSize', 20)
xlabel('Time, s')
ylabel('Normalized amplitude')
title('The signal in the time domain')
% plotting of the spectrum (1) (spectrogram)
figure(2)
%specgram(x, 512, fs, [], 500) % 97% overlapping
spectrogram(x, 512, 500, [], fs, 'yaxis') % 97% overlapping
h = colorbar;
set(h, 'FontName', 'Times New Roman', 'FontSize', 20)
ylabel(h, 'Magnitude, dB');
set(gca, 'FontName', 'Times New Roman', 'FontSize', 20)
xlabel('Time, s')
ylabel('Frequency, Hz')
title('Spectrogram of the signal')
% FFT
win = hanning(N);
K = sum(win)/N; % coherent amplification of the window
X = abs(fft(x.*win)); % fast fourier transform
Xm = X(1:N/2); % getting a first half of the spectrum without Nyquist frequency at N/2+1
Xm = Xm/(N/2); % computing of the amplitudes
Xm(1,1) = Xm(1,1)/2; % correction of the DC component
Xm = Xm/K; % correction due to coherent amplification
% plotting of the spectrum
f = (0:N/2-1)*fs/N;
figure(3)
semilogx(f/1000, 20*log10(Xm))
xlim([min(f/1000) max(f/1000)])
grid on
set(gca, 'FontName', 'Times New Roman', 'FontSize', 20)
title('Amplitude Spectrum of the signal')
xlabel('Frequency, kHz')
ylabel('Amplitude, dBV')
% plotting of the histogram
figure(4)
histfit(x)
xlim([-1.1*max(abs(x)) 1.1*max(abs(x))])
grid on
set(gca, 'FontName', 'Times New Roman', 'FontSize', 20)
xlabel('Signal value')
ylabel('Number of samples')
title('Probability distribution')
set(gca, 'FontName', 'Times New Roman', 'FontSize', 16)
legend('probability distribution of the signal', 'standard normal distribution')
% Sigma and Mu estimation
[u s] = normfit(x);
sigma = num2str(s);
mu = num2str(u);
disp(['Sigma = ' sigma])
disp(['Mu = ' mu])
% computing of the peak (crest) factor
rms = sqrt(mean(x.^2));
x = abs(x);
peak = max(x);
Q = 20*log10(peak/rms);
Qstr = num2str(Q);
disp(['Peak (crest) factor Q = ' Qstr ' dB'])
% computing of the dynamic range
maxval = peak;
x(x==0)=[];
minval = min(x);
D = 20*log10(maxval/minval);
Dstr = num2str(D);
disp(['Dynamic range D = ' Dstr ' dB'])
