Matlab code for comparing two audio (microphone) files

Contents

• Matlab Code:
• signal into left and right, stereo records
• time for right signal only (stereo) for original signals
• find new length
• find new time
• find new frequency
• perform fft and plot of the cut orignal signal
• Zoom in to signal baseline noise signal
• plot two signals
• plot fft compare
• calculate and compare two signals
• corlation

Matlab Code:

Author: Minh Anh Nguyen (minhanhnguyen@q.com) This script will perform the following tasks:

%1. read in two "wav" files, plot the signals of these files
%2. fft
% 3. compare input signals, and display a message on the screen whether or not the signals are the same.
% the MSE must be 0, for both signals are the same.
% The script will read in a stereo data (right channel)

close all;% Close all figures (except those of imtool.)
clear all; % Erase all existing variables.
clc;% Clear the command window.
fontSize = 20;
fontSize1 = 14;
fontSize2 = 12;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Please edit files
% set folder
folder = 'C:\Users\Minh anh\Desktop\microphone-project-demo';

%%Open WAV file:
sound1=audioread('.\good 2.wav');
[y1,fs1]=audioread('.\good 2.wav');

% sound2=audioread('.\good 3.wav');
% [y2,fs2]=audioread('.\good 3.wav');

sound2=audioread('.\upper_channel 2.wav');
[y2,fs2]=audioread('.\upper_channel 2.wav');

% sound2=audioread('.\upper_channel_bide 1.wav');
% [y2,fs2]=audioread('.\upper_channel_bide 1.wav');

% sound2=audioread('.\lower_channel 1.wav');
% [y2,fs2]=audioread('.\lower_channel 1.wav');

% sound2=audioread('.\pin1.wav');
% [y2,fs2]=audioread('.\pin1.wav');

sound1_baseline_noise=audioread('.\good 2_noise.wav');
[y1_basenoise,fs1_basenoise]=audioread('.\good 2_noise.wav');
%
% sound2_baseline_noise=audioread('.\good 3_noise.wav');
% [y2_basenoise,fs2_basenoise]=audioread('.\good 3_noise.wav');

sound2_baseline_noise=audioread('\upper_channel 2_noise.wav');
[y2_basenoise,fs2_basenoise]=audioread('.\upper_channel 2_noise.wav');

% sound2_baseline_noise=audioread('\upper_channel_bide 1_noise.wav');
% [y2_basenoise,fs2_basenoise]=audioread('.\upper_channel_bide 1_noise.wav');

% sound2_baseline_noise=audioread('\lower_channel 1_noise.wav');
% [y2_basenoise,fs2_basenoise]=audioread('.\lower_channel 1_noise.wav');

% sound2_baseline_noise=audioread('.\pin1_noise.wav');
% [y2_basenoise,fs2_basenoise]=audioread('.\pin1_noise.wav');

signal into left and right, stereo records

size(y1);
left1=y1(:,1);
right1=y1(:,2);
[MR1,NR1] = size(right1);
size(y2);
left2=y2(:,1);
right2=y2(:,2);
[MR2,NR2] = size(right2);

time for right signal only (stereo) for original signals

Nl = length(right1); % number of points to analyze
lsl = size (right1); % find the length of the data per second
nbits = 2^(length(right1));
t1l = (0:1:length(right1)-1)/fs1;
fx = fs1*(0:Nl/2-1)/Nl;   %Prepare freq data for plot
T1 = 1/fs1; % period between each sample
N = Nl;

N2 = length(right2);
ls2 = size (right2); % find the length of the data per second
nbits2 = 2^(length(right2));
t2 = (0:1:length(right2)-1)/fs2;
fx2 = fs2*(0:N2/2-1)/N2;   %Prepare freq data for plot
T2 = 1/fs2; % period between each sample

% baseline noise is removed from signals
size(y1_basenoise);
left1_basenoise=y1_basenoise(:,1);
right1_basenoise=y1_basenoise(:,2);
[MR1_basenoise,NR1_basenoise] = size(right1_basenoise);
size(y2_basenoise);
left2_basenoise=y2_basenoise(:,1);
right2_basenoise=y2_basenoise(:,2);
[MR2_basenoise,NR2_basenoise] = size(right2_basenoise);

t1_basenoise = (0:1:length(right1_basenoise)-1)/fs1_basenoise;
t2_basenoise = (0:1:length(right2_basenoise)-1)/fs2_basenoise;


%%cut signal to the same length for comparing
voice1 = right1(fs1*1 : fs1*40);
voice2 = right2(fs2*1 : fs2*40);

find new length

N2x = length(voice2);
N1x = length(voice1);

find new time

t1x = (0:1:length(voice1)-1)/fs1;
t2x = (0:1:length(voice2)-1)/fs2;

find new frequency

f2x = fs2*(0:N2x/2-1)/N2x;
f1x = fs1*(0:N1x/2-1)/N1x;

perform fft and plot of the cut orignal signal

NFFT1x = 2 ^ nextpow2(N1x);
Y1x = fft(voice1, NFFT1x)/ N1x;
f1xx = (fs1/ 2 * linspace(0, 1, NFFT1x / 2+1))'; % Vector containing frequencies in Hz
STFFT1x = ( 2 * abs(Y1x(1: NFFT1x / 2+1))); % Vector containing corresponding amplitudes

NFFT2x = 2 ^ nextpow2(N2x);
Y2x = fft(voice2, NFFT2x) / N2x;
f2xx = (fs2 / 2 * linspace(0, 1, NFFT2x / 2+1))'; % Vector containing frequencies in Hz
STFFT2x = ( 2 * abs(Y2x(1: NFFT2x / 2+1))); % Vector containing corresponding amplitudes

% plot of the cut orginal signal
figure; subplot (2,2,1); plot(t1x, voice1);xlim([0,40]); ylim([-2,2]);
%str1=sprintf('Correct loading  with sampling rate = %d Hz, number sample = %d', fs1, N1x);
str1=sprintf('Correct loading signal');
title(str1,'Fontsize', fontSize2);
xlabel('Time (sec)','Fontsize', fontSize1); ylabel('Amplitude','Fontsize', fontSize1); grid on;
subplot (2,2,2); plot(t2x, voice2); xlim([0,40]); ylim([-2,2]);
%str2=sprintf(' User loading with sampling rate = %d Hz, number sample = %d', fs2, N2x);
str2=sprintf('User loading signal');
title(str2,'Fontsize', fontSize2); xlabel('Time (sec)','Fontsize', fontSize1); ylabel('Amplitude','Fontsize', fontSize1); grid on;

% plot of the fft cut signal2
%subplot (2,2,3); plot(f1xx, STFFT1x); xlim([0,1500]); ylim([0,.09]);
figure;
plot(f1xx, STFFT1x); xlim([0,1500]); ylim([0,.09]);
str1ft1=sprintf('Single-sided Power spectrum for Correct loading');
title(str1ft1,'Fontsize', fontSize2);
ylabel ('Magnitude','Fontsize', fontSize1); xlabel ('Frequency [Hz]','Fontsize', fontSize1); grid on;

% display markers at specific data points on DFT graph
indexmin1 = find(min(STFFT1x) == STFFT1x);
xmin1 = f1xx(indexmin1);
ymin1 = STFFT1x(indexmin1);
indexmax1 = find(max(STFFT1x) == STFFT1x);
xmax1 = f1xx(indexmax1);
ymax1 = STFFT1x(indexmax1);
strmax1 = ['Frequency = ',num2str(xmax1),' Hz','    ','Magnitude=', num2str(ymax1),''];
text(xmax1,ymax1,strmax1,'HorizontalAlignment','Left');
%subplot (2,2,4); plot(f2xx, STFFT2x); xlim([0,1500]); ylim([0,.09]);
figure; plot(f2xx, STFFT2x); xlim([0,1500]); ylim([0,.09]);
str1ft2=sprintf('Single-sided Power spectrum for User loading');
title(str1ft2,'Fontsize', fontSize2);
ylabel ('Magnitude', 'Fontsize', fontSize1); xlabel ('Frequency [Hz]','Fontsize', fontSize1); grid on;
% display markers at specific data points on DFT graph
indexmin2 = find(min(STFFT2x) == STFFT2x);
xmin2 = f2xx(indexmin2);
ymin2 = STFFT2x(indexmin2);

indexmax2 = find(max(STFFT2x) == STFFT2x);
xmax2 = f2xx(indexmax2);
ymax2 = STFFT2x(indexmax2);
strmax2 = ['Frequency = ',num2str(xmax2),' Hz','    ','Magnitude =', num2str(ymax2),' '];
text(xmax2,ymax2,strmax2,'HorizontalAlignment','Left');

  

 

Zoom in to signal baseline noise signal

figure; subplot(3,2,1); plot(t1_basenoise,right1_basenoise); xlim([0,10]); ylim([-2,2]);
str1=sprintf('Correct loading with Zoom in');
title(str1,'Fontsize', fontSize1);
xlabel('Time (secs)','Fontsize', fontSize1); ylabel('Amplitude','Fontsize', fontSize1); grid on;

subplot(3,2,3); plot(t2_basenoise,right2_basenoise,'r'); xlim([0,10]); ylim([-2,2]);
str2=sprintf('User loading with Zoom in');
title(str2, 'Fontsize', fontSize1);
xlabel('Time (secs)','Fontsize', fontSize1); ylabel('Amplitude','Fontsize', fontSize1); grid on;

 

plot two signals

figure, plot(t1_basenoise,right1_basenoise, t2_basenoise,right2_basenoise,'r'); xlim([0,10]); ylim([-2,2]);
str2=sprintf('Comparing two loading signals with Zoom in ');
title(str2, 'Fontsize', fontSize1);
xlabel('Time (secs)','Fontsize', fontSize1); ylabel('Amplitude','Fontsize', fontSize1); grid on;
legend ('Correct loading signal', 'User loading signal', 'Fontsize', fontSize);

Warning: Using an integer to specify the legend location is not supported.
Specify the legend location with respect to the axes using the 'Location'
parameter. 
Warning: Ignoring extra legend entries. 

plot fft compare

figure, plot(f1xx, STFFT1x,f2xx, STFFT2x,'r'); xlim([0,1500]); ylim([0,.09]);
str1ft3=sprintf('Comparing FFT for two loading signals');
title(str1ft3,'Fontsize', fontSize1);
ylabel ('Magnitude', 'Fontsize', fontSize1); xlabel ('Frequency [Hz]','Fontsize', fontSize1); grid on;

indexmax3 = find(max(STFFT2x) == STFFT2x);
xmax3 = f2xx(indexmax3);
ymax3 = STFFT2x(indexmax3);
strmax3 = ['Frequency = ',num2str(xmax3),' Hz','    ','Magnitude= ', num2str(ymax3),' '];
legend ('Correct loading signal', 'User loading signal','Fontsize', fontSize);

Warning: Using an integer to specify the legend location is not supported.
Specify the legend location with respect to the axes using the 'Location'
parameter. 
Warning: Ignoring extra legend entries. 

calculate and compare two signals

Calculate the MSE

D= voice1 - voice2;
MSE=mean(D.^2);
    %msgbox(strcat('MSE value is= ',mat2str(),' MSE'));
    %msgbox(strcat('MSE value is= ',mat2str(MSE)));
%if MSE ==0;
 %if MSE <= 0.00135;  %%v this will work
     if MSE <= 5.2731e-09 ;  %%v this will work
    disp('Signal are the same');
else
%        msgbox(strcat('MSE value is= ',mat2str(MSE)));
         disp('Signal are not the same');
 msgbox('Centrifuge loading error. Please check bearings are seated properly');
     end


 % %% play back play back the audio files
% disp('Playing (correctly loaded).');
% sound(y1,fs1);
% pause; %create a 2 second pause
%
%
% disp('Playing the user loading.');
% sound(y2,fs2);
% pause





%%Calculate variance, mean and standard deviation of sounds

D1c= STFFT1x - STFFT2x;
MSEfft=mean(D1c.^2);
disp(['MSEfft' num2str(MSEfft)]);

mean1= sum(right1)/N;
mean2= sum(right2)/N2;
stdev1 = sqrt ((sum(right1-mean1).^2)/2);
stdev2 = sqrt ((sum(right2-mean2).^2)/2);



% var(right1(1:length(right2))-right2) %for the whole signal y1-y2
var(voice1(1:length(voice2))-voice2); %for the voice

 fprintf('digital statitics of voice 1 a same length signal\n\n');
 fprintf('mean of voice 1:%f\n', mean(voice1));
 fprintf('standard deviation of voice 1:%f\n', std(voice1));
 fprintf('variance of voice 1:%f\n', std(voice1)^2);
 fprintf('average power of voice 1:%f\n', mean(voice1.^2));
 fprintf('average magnitude of voice 1:%f\n', mean(abs(voice1)));
%  prod1 = sound1(1:N-1) * sound1(2:N);
%  crossing1 = length (find(prod1<0));
%  fprintf(' zero crossing of sound 1:% 0.f\n', crossing1);
%
 fprintf('\ndigital statitics of voice 2\n\n');
 fprintf('mean of voice 2:%f\n', mean(voice2));
 fprintf('standard deviation of voice 2:%f\n', std(voice2));
 fprintf('variance of voice 2:%f\n', std(voice2)^2);
 fprintf('average power of voice 2:%f\n', mean(voice2.^2));
 fprintf('average magnitude of voice 2:%f\n', mean(abs(voice2)));
%  prod1 = sound1(1:N-1) * sound1(2:N);
%  crossing1 = length (find(prod1<0));
%  fprintf(' zero crossing of sound 1:% 0.f\n', crossing1);
%
% %  var(right1(1:length(right2))-right2) %for the whole signal y1-y2
%  var(voice1(1:length(voice2))-voice2); %for the voice
%
%

Signal are not the same
MSEfft5.7938e-08
digital statitics of voice 1 a same length signal

mean of voice 1:-0.000105
standard deviation of voice 1:0.251742
variance of voice 1:0.063374
average power of voice 1:0.063374
average magnitude of voice 1:0.197595

digital statitics of voice 2

mean of voice 2:0.000030
standard deviation of voice 2:0.235003
variance of voice 2:0.055227
average power of voice 2:0.055227
average magnitude of voice 2:0.185665

correlation

[C1, lag1] = xcorr(abs((fft(voice1))),abs((fft(voice2)) ));
 figure, plot(lag1/fs1,C1);
    ylabel('Amplitude'); grid on
    title('Cross-correlation between Correct and User  loading files')

% %% plot orginal signal. for debug. %% team does not like long length
% %% perform fft and plot of orignal signal
% %
%  figure; subplot (2,2,1); plot(t1, sound1);
%  title ('Signal of sound 1 , frequency = 44100 Hz and number sample = 2742270');
%  xlabel('time (sec)'); ylabel('relative signal strength'); grid on;
% %
%  subplot (2,2,2); plot(t2, sound2);
%   title ('Signal of sound 2, frequency = 44100 Hz and number sample = 2674665');
%  xlabel('time (sec)'); ylabel('relative signal strength'); grid on;


% %% fft of orignal signal
% NFFT1 = 2 ^ nextpow2(N);
% Y1 = fft(y1, NFFT1)/ N;
% f1 = (fs1 / 2 * linspace(0, 1, NFFT1 / 2+1))'; % Vector containing frequencies in Hz
% STFFT1 = ( 2 * abs(Y1(1: NFFT1 / 2+1))); % Vector containing corresponding amplitudes
% subplot (2,2,3); plot(f1, STFFT1); title ('Single-sided Power spectrum for Sound1');
% %subplot (2,2,3); plot(abs(fft(sound1))); title ('plot fft of Sound 1'); %plot Magnitude of fft
% ylabel ('Magnitude |y(f)|'); xlabel ('frequency [Hz]'); grid on;
%
% NFFT2 = 2 ^ nextpow2(N2);
% Y2 = fft(y2, NFFT2) / N2;
% f2 = (fs2 / 2 * linspace(0, 1, NFFT2 / 2+1))'; % Vector containing frequencies in Hz
% STFFT2 = ( 2 * abs(Y2(1: NFFT2 / 2+1))); % Vector containing corresponding amplitudes
% subplot (2,2,4); plot(f2, STFFT2); title ('Single-sided Power spectrum for Sound2');
% % subplot (2,2,4); plot(abs(fft(sound2))); title ('plot fft of Sound2');
% ylabel ('Magnitude |y(f)|'); xlabel ('frequency [Hz]'); grid on;

Published with MATLAB® R2016a