Reading & Plotting Audio Files in MATLAB

Vladimir Kulyukin

---

Problem
  

This post describes how to read in an audio file and plotting it in MATLAB. We will work with this audio file (calcium.wav) that has a recording of the word calcium. We will assume, for the sake of simplicity, that the file is saved in C:\Users\Vladimir\audio_files. 

---

Reading Audio Files and Computing its Statistics

Let us add the path to the MATLAB's path list to reduce the amount of typing, read the file in, and play it.

addpath C:\Users\Vladimir\audio_files
[speech_vector, freq, bits_per_sample] = wavread('calcium.wav');
soundsc(speech_vector, freq);

The variable speech_vector is a digital speech vector. We can compute its length.

>> length(speech_vector)
ans =
      124928

There are 124928 samples in the vector. The second variable freq is the frequency. Recall that freq = 1/T where T is the time duration of the sample. Let us find out what the value of freq is.

>> freq  
freq =
       44100

So the frequency is 44100 Hz. The third variable, bits_per_sample, tells us the number of bits per sample. In our case, it is equal to 16.

>> bits_per_sample
bits_per_sample =
    16

Let us compute the duration of each sample T in milliseconds using the formula T = 1/freq. 

>> 1.0/freq
ans =
   2.2676e-05
 
Let us compute the duration of each sample in seconds.

>> 1.0/freq*1000
ans =
    0.0227 

Here is how we can compute the duration of the entire file in seconds.

>> length(speech_vector)*1.0/freq*1000
ans =
   2.8328e+03 

So, the duration of the entire file is 2.833 seconds. Figure 1 shows the file plotted in Audacity. The time line above confirms that the duration of the file is 2.83 seconds.

Figure 1. Calcium.wav plotted in Audacity

A shorter way of computing the time duration of the file is to divide the length by the speech vector by its frequency.

>> length(speech_vector)/freq
ans =
    2.8328

We can play the file faster or slower. For example, this command plays the audio file twice as fast.

>> soundsc(speech_vector, 2*freq)

This command plays the audio file twice as slow.

>> soundsc(speech_vector, freq/2.0)

---

Plotting Audio Files

Figure 2 shows how we can plot the speech vector read with the audio file with the command plot(speech_vector).

Figure 2. Output of the command plot(speech_vector)

We can plot the speech vector to bring it in line with the Audacity plot in Figure 1 by changing in the x-axis to seconds. Here is how.

step = 1/length(speech_vector);
xval = xval*(length(speech_vector)/freq);
plot(xval, speech_vector);
 
Figure 3 gives the output of the above instructions.

Figure 3. Changing the x-axis to seconds

The same plot can be produced with the following command, which is not as clear.

 plot((0:1/length(speech_vector):1-1/length(speech_vector))*(length(speech_vector)/freq), speech_vector);

Here is how we can label the x- and y-axis, add a grid, a legend, and save the plot as a png file.

plot(xval, speech_vector);
xlabel('Time [seconds]');
ylabel('Amplitude');
grid on;
legend('Calcium.wav');
figure(1), print -dpng 'C:\Users\Vladimir\audio_files\Calcium.png'; 

---

Absolute Magnitude of the FFT of Sinusoids with Unit Amplitudes and Unit Periods

Vladimir Kulyukin

---

Problem

This post briefly describes how to obtain and plot the absolute magnitude of a sinusoid with a unity amplitude and period. and normalized frequencies of the FFT. Figure 1 sketches a proof of why the absolute magnitude of the FFT of length N of  a sinusoid with a unit amplitude and period  is equal to N/2.
  

Figure 12. Plots of sinusoid 04, its absolute magnitudes, and normalized frequencies

Figure 12. Plots of sinusoid 04, its absolute magnitudes, and normalized frequencies

Figure 1. Absolute Magnitude of the FFT

---

Sinusoid 01
 

Figure 2.  Sinusoid 01 with A=1, T=1, L=4, x=1

Here is the MATLAB code that encodes sinusoid 01 defined in Figure 2.

N=64;       %% this should be a power of 2
T=1;        %% this is the period
A=1;        %% amplitude
phi=0.0;    %% phase of the sinusoid
fs = 1/T;   %% fs
t = 0:N-1   %% discrete time axis
L = 4;
x = 1;
wx = x*2*pi*fs/L %% wx frequency
s_01a = A*cos(wx*t + phi); %% input signal
FFT_01a = fft(s_01a); %% spectrum
%% Plot the sinusoid
figure(1);
subplot(3,1,1);
plot(t, s_01a, '*k');
ti = 0:.1:N-1; %% interpolate time axis in 0.1 increments
hold on;
plot(ti, A*cos(wx*ti + phi), '-k'); grid off;
title('Sinusoid signal: Example 01a');
xlabel('Time (samples)');
ylabel('Amplitude');
text(-8, 1, 'a)');
hold off;
%% Plot absolute magnitude at normalized frequencies
magFFT_01a = abs(FFT_01a);
disp(magFFT_01a);
tn = 0:1/N:1-1/N;
subplot(3, 1, 2);
stem(tn, magFFT_01a, 'ok'); grid on;
xlabel('Frequency (cycles per sample)');
ylabel('Magnitude (linear)');
text(-.11, 40, 'b');
%% Plotting dB scales
dbSpectrum_01a = 20*log10(magFFT_01a); %% convert into dB
subplot(3, 1, 3);
plot(tn, dbSpectrum_01a, '--ok'); grid on;
axis([0 1 -350 50]);
xlabel('Normalized frequency (cycles per sample)');
ylabel('Magnitude (dB)');
text(-.15, 50, 'c)');


Figure 3 gives the graphs computed and plotted by the above MATLAB code.

Figure 3. Plots of sinusoid 01, its absolute magnitudes, and normalized frequencies

---

Sinusoid 02

Figure 7. Sinusoid 02 with A=1, T=1, L=8, x=1

Here is the MATLAB code for sinusoid 02 defined in Figure 7.

N=64;       %% this should be a power of 2
T=1;        %% this is the period
A=1;        %% amplitude
phi=0.0;    %% phase of the sinusoid
fs = 1/T;   %% fs
t = 0:N-1   %% discrete time axis
L = 8;
x = 1;
wx = x*2*pi*fs/L %% wx frequency
s_01b = A*cos(wx*t + phi); %% input signal
FFT_01b = fft(s_01b); %% spectrum
%% Plot the sinusoid
figure(1);
subplot(3,1,1);
plot(t, s_01b, '*k');
ti = 0:.1:N-1; %% interpolate time axis in 0.1 increments
hold on;
plot(ti, A*cos(wx*ti + phi), '-k'); grid off;
title('Sinusoid signal: Example 02');
xlabel('Time (samples)');
ylabel('Amplitude');
text(-8, 1, 'a)');
hold off;
%% Plot absolute magnitude at normalized frequencies
magFFT_01b = abs(FFT_01b);
disp(magFFT_01b);
tn = 0:1/N:1-1/N;
subplot(3, 1, 2);
stem(tn, magFFT_01b, 'ok'); grid on;
xlabel('Frequency (cycles per sample)');
ylabel('Magnitude (linear)');
text(-.11, 40, 'b');
%% Plotting dB scales
dbSpectrum_01b = 20*log10(magFFT_01b); %% convert into dB
subplot(3, 1, 3);
plot(tn, dbSpectrum_01b, '--ok'); grid on;
axis([0 1 -350 50]);
xlabel('Normalized frequency (cycles per sample)');
ylabel('Magnitude (dB)');
text(-.15, 50, 'c)');

Figure 8 gives the graphs computed and plotted by the above MATLAB code.

Figure 8. Plots of sinusoid 02, its absolute magnitudes, and normalized frequencies

---

Sinusoid 03
 

Figure 9. Sinusoid with A=1, T=1, L=16, x=1

Here is the MATLAB code for sinusoid 03 in Figure 9.

N=64;       %% this should be a power of 2
T=1;        %% this is the period
A=1;        %% amplitude
phi=0.0;    %% phase of the sinusoid
fs = 1/T;   %% fs
t = 0:N-1   %% discrete time axis
L = 16;
x = 1;
wx = x*2*pi*fs/L %% wx frequency
s_01c = A*cos(wx*t + phi); %% input signal
FFT_01c = fft(s_01c); %% spectrum
%% Plot the sinusoid
figure(1);
subplot(3,1,1);
plot(t, s_01c, '*k');
ti = 0:.1:N-1; %% interpolate time axis in 0.1 increments
hold on;
plot(ti, A*cos(wx*ti + phi), '-k'); grid off;
title('Sinusoid signal: Example 03');
xlabel('Time (samples)');
ylabel('Amplitude');
text(-8, 1, 'a)');
hold off;
%% Plot absolute magnitude at normalized frequencies
magFFT_01c = abs(FFT_01c);
disp(magFFT_01c);
tn = 0:1/N:1-1/N;
subplot(3, 1, 2);
stem(tn, magFFT_01c, 'ok'); grid on;
xlabel('Frequency (cycles per sample)');
ylabel('Magnitude (linear)');
text(-.11, 40, 'b');
%% Plotting dB scales
dbSpectrum_01c = 20*log10(magFFT_01c); %% convert into dB
subplot(3, 1, 3);
plot(tn, dbSpectrum_01c, '--ok'); grid on;
axis([0 1 -350 50]);
xlabel('Normalized frequency (cycles per sample)');
ylabel('Magnitude (dB)');
text(-.15, 50, 'c)');

Figure 10 gives the graphs computed and plotted by the above MATLAB code.
 

Figure 10. Plots of sinusoid 03, its absolute magnitudes, and normalized frequencies

---

Sinusoid 04
 

Figure 11. Sinusoid 04 with A=1, T=1, L=4, x=2

Here is the MATLAB source for  sinusoid 04 defined in Figure 11.

N=64;       %% this should be a power of 2
T=1;        %% this is the period
A=1;        %% amplitude
phi=0.0;    %% phase of the sinusoid
fs = 1/T;   %% fs
t = 0:N-1   %% discrete time axis
L = 4;
x = 2;
wx = x*2*pi*fs/L %% wx frequency
s_01d = A*cos(wx*t+phi); %% input signal
FFT_01d = fft(s_01d); %% spectrum
%% Plot the sinusoid
figure(1);
subplot(3,1,1);
plot(t, s_01d, '*k');
ti = 0:.1:N-1; %% interpolate time axis in 0.1 increments
hold on;
plot(ti, A*cos(wx*ti + phi), '-k'); grid off;
title('Sinusoid input: Example 04');
xlabel('Time (samples)');
ylabel('Amplitude');
text(-8, 1, 'a)');
hold off;
%% Plot absolute magnitude at normalized frequencies
magFFT_01d = abs(FFT_01d);
disp(magFFT_01d);
tn = 0:1/N:1-1/N;
subplot(3, 1, 2);
stem(tn, magFFT_01d, 'ok'); grid on;
xlabel('Frequency (cycles per sample)');
ylabel('Magnitude (linear)');
text(-.11, 40, 'b');
%% Plotting dB scales
dbSpectrum_01d = 20*log10(magFFT_01d); %% convert into dB
subplot(3, 1, 3);
plot(tn, dbSpectrum_01d, '--ok'); grid on;
axis([0 1 -350 50]);
xlabel('Normalized frequency (cycles per sample)');
ylabel('Magnitude (dB)');
text(-.15, 50, 'c)');

Figure 12 gives the graphs computed and plotted by the above MATLAB code.

Figure 12. Plots of sinusoid 04, its absolute magnitudes, and normalized frequencies

---