How do I go about intentionally aliasing a signal, and upsampling work with audio file in MATLAB

Question

My objective is to show the results of:

• How do I go about intentionally aliasing a signal?
• undersampling
• Nyquist rate sampling and oversampling

I'm first going to get the MATLAB code working well, then I'll rewrite it for the MSP430.

First I need to downsample the .wav file to get an incomplete or impartial data stream that I can then reconstuct. Here the flow chart:
_{analog -> sampling analog filter -> ADC -> resample down -> resample up -> DAC -> reconstruction analog filter}

This is my first signal processing project using MATLAB (v2010), so be nice. Original code below; updated code here.

    %Play decimated file ( soundsc(y,fs) ) 
%Play Original file ( soundsc(play,fs ) )
%Play reconstucted File ( soundsc(final,fs) )

[piano,fs]=wavread('piano.wav'); % loads piano
play=piano(:,1); % Renames the file as "play"

t = linspace(0,time,length(play));          % Time vector
x = play;
y = decimate(x,25);

stem(x(1:30)), axis([0 30 -2 2])   % Original signal
title('Original Signal')
figure
stem(y(1:30))                        % Decimated signal
title('Decimated Signal')

%changes the sampling rate

fs1 = fs/2;
fs2 = fs/3;
fs3 = fs/4;
fs4 = fs*2;
fs5 = fs*3;
fs6 = fs*4;

wavwrite(y,fs/25,'PianoDecimation');


%------------------------------------------------------------------

%Downsampled version of piano is now upsampled to the original
[PianoDecimation,fs]=wavread('PianoDecimation.wav'); % loads piano
play2=PianoDecimation(:,1); % Renames the file as "play

%upsampling
UpSampleRatio = 2;  % 2*fs = nyquist rate sampling
play2Up=zeros(length(PianoDecimation)*UpSampleRatio, 1);
play2Up(1:UpSampleRatio:end) = play2; % fill in every N'th sample

%low pass filter

ResampFilt = firpm(44, [0 0.39625 0.60938 1], [1 1 0 0]);


fsUp = (fs*UpSampleRatio)*1;
wavwrite(play2Up,fsUp,'PianoUpsampled');

%Plot2
%data vs time plot
time=(1/44100)*length(play2);
t=linspace(0,time,length(play2));
stem(t,play2)
title('Upsampled graph of piano')
xlabel('time(sec)');
ylabel('relative signal strength')



[PianoUpsampled,fs]=wavread('PianoUpsampled.wav'); % loads piano
final=PianoUpsampled(:,1); % Renames the file as "play"


%-------------------------------------------------------------
%resampleing
[piano,fs]=wavread('piano.wav'); % loads piano
x=piano(:,1); % Renames the file as "play"
m = resample(x,3,2);

**%this is were i would need to sample it at the different frequecys (both above and below and at) nyquist frequency.*I think.***

soundsc(left,fs) % shows the resaultant audio file , which is the same as original ( only at or above nyquist frequency however) 

How can I improve this? How do I sample at different frequencies?

^{Note that I've crossposted this on SO, here.}

Answer 1

If you want to downsample in a way that produces aliasing, just throw away some fraction of the samples, ie, only keep every 2nd sample, etc.

If you want to downsample in a way that doesn't alias, precede the decimation with a lowpass filter having a cutoff below the nyquist frequency of the sample rate you'll have after you decimate. Since decimating involves throwing away samples and an FIR lowpass filter has no future dependence on its results, you don't have to bother doing a convolution dot product for the outputs you'd only discard in the decimation, instead you can just skip to the next one.

To upsample, insert zeroes between the input samples (use zeroes, do not repeate the input values). This produces alias frequencies, so follow with a lowpass filter below the nyquist rate of the original sample rate. Again, you can design the convolution computation to skip over the zero inputs.

For small integer ratios of sample rates, discarding or zero inserting samples and FIR filters work very well. For larger ratios, CIC filters may be preferred as they can be efficiently computed in hardware using wrapping binary computation (you intentionally let the calculations overflow!), though they have a sinc function-shaped rather than rectangular frequency response and so need a clean up filter to select only the rectangular area near the tip of the sinc response and/or compensate for the droop of the sinc.

Answer 2

I am going to start simple, and build up.

F= Frequency

F(Hz=1/s) E.x. 100Hz = 1000 (Cyc/sec) F(s)= 1/(2f)

Example problem: 1000 hz = Highest frequency 1/2(1000hz) = 1/2000 = 
5x10(-3) sec/cyc or a sampling rate of 5ms

If your highest frequency is 1000Hz and you sample with a period of 5mS(500Hz) your signal has already lost a large amount of its information. You are sampling at 1/4th of the nyquist rate. For a 1000Hz wave, you must sample at at least 2000Hz, and that is in a perfect world. You really need to go a bit faster than this due to technology limitations, but this is beyond what you need to know.

What is my nyquist frequency?

Wait, do you mean Nyquist Frequency or Nyquist rate?

This is probably seeming a pretty trivial distinction, but many people are referencing different things when they switch terms, some people can control the "Symbol Rate" some can control the "Sampling Rate," some can control both(lucky engineers).

You need to sample at a frequency that is greater than double your highest frequency.

If you go lower than this your signal will alias, this is irreversible once it has happened, the best you can do is filter off the corrupted part of your signal, but this does not recover information, just removes the corrupted frequencies.

Downsampling in Matlab

Lets say that you have your signal, you have avoided aliasing, and you want to downsample to half the points. This is very easy in matlab. Take a look at the decimate function. This will handle the mess of making sure your down-sampling does not introduce aliasing, at a cost. When it filters off the higher frequencies, they are lost. This is what you ask for when you half your datapoints. You can just use the downsample function. If you had the higher frequencies, it will alias. If you know there should not be those frequencies there, you can always use decimate to make sure that noise does not sneak into your system.

I would like to go into more detail, but I had to cut myself short as I am out of time for today. Leave me a comment if there is more that would help you, and I will come back and add more information. I am sure I did not cover everything you wanted.