Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
144 lines (101 sloc) 3.52 KB
close all
clc
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
% ALIASING %
% Author: Cecilia Casarini %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% In this script the concept of aliasing is illustrated by plotting a
% signal sampled above and below the Nyquist frequency and by playing the
% resulting signals. The script is interrupted by pauses in order to have
% more time to think about every step. Click any key to continue.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Building the "analog" signal (it is not really analog, of course, but it
% has many points, so it resembles a continuous signal
% Number of points
NP = 50000;
% Signal length (s) -increase tmax if you want to hear a longer signal(but
% you will have a more dense graph.)
tmin = 0;
tmax = 0.01;
% Time axis
tc = linspace(tmin, tmax, NP);
% Frequency of the signal (Hz)
f = 1000;
% Signal
sc = sin(2*pi*f*tc);
% Plot the signal s
plot(tc, sc)
title('Continuous Signal');
xlabel('Time (s)');
ylabel('Amplitude');
hold on
pause
% 2. Sampling the signal at 44100 Hz, which is more than 2 times the
% frequency. Therefore the signal will have enough points to represent the
% correct frequency.
% Sample rate (Hz)
Fs1 = 44100;
% Sampling period
Ts1 = 1/Fs1;
% Time line
t1 = tmin:Ts1:tmax;
% Signal
s1 = sin(2*pi*f*t1);
% Plotting the samples
plot(t1,s1, '.')
title('Correctly Sampled Signal');
xlabel('Time (s)');
ylabel('Amplitude');
hold off
pause
% Playing the resulting frequency
soundsc(s1, Fs1)
% Saving the .wav file
audiowrite('Original_1000Hz.wav',s1,Fs1)
pause
% 3. Sampling the signal at 1500 Hz, which is less than 2 times the
% frequency. Therefore the signal will be aliased. The frequency heard will
% be the difference between the sampling rate and the original frequency.
% In this case 1500 Hz - 1000 Hz = 500 Hz, which is an octave below the
% original frequency. It can be heard clearly (if you know how an octave
% interval sounds like, or also easily visualised: when the 500 Hz signal
% will have completed one cycle, the 1000 Hz will have completed 2 cycles.
% Sample rate (Hz)
Fs2 = 1500;
% Sampling period
Ts2 = 1/Fs2;
% Time line
t2 = tmin:Ts2:tmax;
% Signal
s2 = sin(2*pi*f*t2);
% Plotting again the continuous signal
plot(tc, sc)
hold on
% Plotting the sampled signal at Fs = 1500 Hz
plot(t2,s2, '.')
title('Aliased Signal');
xlabel('Time (s)');
ylabel('Amplitude');
legend('Continuous Original', 'Sampled Aliased');
pause
% Showing that the aliased frequency is indeed 500 Hz instead
% of 1000 Hz
% Shannon-Nyquist theorem provides a way to calculate the aliased frequency
k = round(-(f/Fs2+1/2));
fa = f+k*Fs2;
% Aliased Signal
s3 = sin(2*pi*(fa)*tc);
% Plot
plot(tc,s3)
legend('Continuous Original', 'Sampled Aliased','Continuous Aliased');
hold off
% Playing the aliased frequency (you can hear it's an octave below!)
pause
soundsc(s2, Fs2)
% Saving the .wav file - whoops, it won't save the following, probably
% matlab has an antialiasing filter?
% audiowrite('Aliased_500Hz.wav',s2,Fs2)
% Well... I'll save s3 instead! We already heard with soundsc that it's
% indeed an octave below
audiowrite('Aliased_500Hz.wav',s3,NP)