Useful Matlab Functions and Scripts for Audio Signals and Systems EE513

This page contains mfile and pfile downloads used in class examples and assignments. Mfiles are simply text(ascii) files. They can be copied and pasted into the Matlab editor or workspace, or downloaded to you computer from the links below. Pfiles are parsed mfiles that cannot be edited or their contents seen. Most mfiles were hastily written by Kevin D. Donohue, so some debugging may be required.

tone.m  This function generates a tone signal at a given frequency and time duration.

tonec.m  This function generate a complex tone (simultaneous discrete frequencies) for a given vector of frequencies and relative weights.

convolution_show.m  This script generates a series of plots to graphically illustrate the discrete convolution between an impulse response and an input signal to compute the output.

freqres.m  This script plots the frequency response and impulse response for a system with a Z-transform given by H(z) = 1 / (1 - .5*z^-1).

aliasexm.m  This script creates a movie to how frequencies greater than half the sampling rate get aliased.

filter_analysis.m  This script shows an example of analyzing a filter/system with the magnitude response, phase response, and group delay.  These are computed directly with basic Matlab operations and also using the Matlab's function freqz and grpdelay for comparison.

fir1ex.m  This script solves the example in lecture notes comparing filter orders for the linear phase FIR filters (designed using the windowing method).  High-pass filters of orders of 10 and 50 are designed with a cut-off of 500 Hz.  The freqz command is used to examine the frequency response (phase and magnitude), the grpdelay command is used to examine the group delay for the filter, and finally a frequency swept signal (from 20Hz  to 2000Hz is generated and filtered with both filters for comparisons.

chebyex.m  This script solves the example in lecture notes comparing filter orders for the Chebyshev type 2 IIR filters.  High-pass filters of orders 5 and 10 are designed with a cut-off of 500 Hz and a stopband ripple of 30 dB down. The freqz command is used to examine the frequency response (phase and magnitude), the grpdelay command is used to examine the group delay for the filter, and finally a frequency swept signal (from 20Hz to 2000Hz is generated and filtered with both filters for comparisons.

ultsigfil.m  This script applies a Butterworth filter to an ultrasonic signal corrupted by out-of-band noise from the instrumentation electronics.  The signal consists of a series of echoes from tissue structures from a 7 MHz probe.  The DFT magnitude of the raw signal is plotted to examine where the energy bands are for the signal and noise. From this figure it is decided that a bandpass filter from 4MHz to 8.5 MHz should be used to pass the signal energy and attenuated the noise. A Butterworth filter of order 5 used and the signal is filtered forwards and backwards to ensure zero phase distortion.  The script requires data in mat file ultdata.mat.

voicefilt.m  This script presents a filtering example from the lecture notes where a sample voice signal is given (random speech located in wave file nad1.wav) and the average spectrum of this sound is used to design a filter to a different voice signal from the same person (located in wave file nad3.wav) corrupted by additive white noise.  The script requires data in wave files: nad1.wav and nad1.wav.

myvoicefilt.m  This script presents a filtering example where a sample voice signal is given (a repeated vowel sound located in wave file myvowel.wav) and the average spectrum of this sound is used to design a filter to separate the voice signal (located in wave file mytalk.wav) from white noise.  The script requires data in wave files: myvowel.wav and mytalk.wav.

proneyex.m - This script solves the lecture note example of deriving an ARMA model of a room impulse response from the data using Proneys  method.  Wave file recording of a clap sound in the room is required to run this script: offroomimp2.wav .

sincex1.m - Sinc function weighting example steps through points of a sampled rectangular pulse to interpolate points that already exist (based on the original sampling). This demonstrates that the sinc interpolation on existing points simply returns the point back again, noticing in the graphics of this script that the sinc nulls fall on all the other existing samples, resulting in a zero contribution to the interpolated point.

sincex2.m - Sinc function weighting example steps through points of a sampled rectangular pulse to interpolate points between existing points (based on a higher sampling frequency). This demonstrates sinc weighting on neighboring points to interpolate points between existing points.  In the graphics of this script notice the intersection of the sinc with all the other existing samples, resulting in a weighted contribution to the interpolated point.

sincinterp.m - Sinc function interpolation example upsamples a rectangular pulse through interpolation.  In the final plot note the effects of the bandlimited nature of the sinc function interpolation (recall it is the FT pair with the ideal low pass filter).  The overshoot and undershoot around the edges is sometimes referred to as Gibbs phenomenon.

rayring.m - This script simulates a tone (ring) with a Rayleigh envelope.  The signal is truncated and the FFT is taken in order to examine its spectral content.  The script can be used to study the effects of zero padding, signal truncation, and signal duration on the estimated spectrum.

spectroex.m - This script simulates a frequency sweep signal and takes its spectrogram to examine its spectral content changes over time, and compares spectrogram result to the FFT of the whole sweep signal.  In addition, parameters and settings such as zero padding, window size window shape, window overlap, as well as signal characteristics can be changed to observe effects.

colorednoise.m - This script generates colored noise by filtering white noise with a Butterworth band-pass filter. Figures are then generated to compare the PSD of the colored noise to the filter transfer function magnitude.  In addition, the autocorrelation is taken of the input white noise sequence and the output colored noise sequence for comparisons.

reverbex.m - This script loads in a recorded clap noise (download from clap.wav) and simulates a reverberation based on 3 fundamental delay intervals.  It shows how to use the Hilbert transform to compute the signal envelope and estimate the RT60 time.  It also applies the autocorrelation function on both the envelope and original reverberated signal to look for peaks corresponding to the fundamental delays of the reverb.  To run this script you will need the function to implement the reverb with multiple delays, mrevera.m, and the function to generate filter coefficients for the all-pass reverb, areverb.m (Optional download: If you comment in and out certain lines in mrevera.m you can apply a plain reverb that has a comb-like magnitude response.  For this you will need the plain reverb function preverb.m.)

ampsqa.p - This function simulates digitizing a signal with nonlinear (amplifier) characteristics, a 16 bit quantization over a ± 3 volt range, and a sampling rate of 16kHz. The input should be a signal sampled at a rate higher than or equal to 16kHz for the simulation to be effective (if a lower rate is used, it will be upsampled, which is not realistic for an analog to digital converter).  The function is called by:

>>    sigout = ampsqa(sigin, fso)

where SIGIN represents the analog signal (not quantized) and created with a sampling rate of FSO Hz. The signal is passed through an anti-aliasing filter and sampled at 16kHz (resamples).  The amplifier has a high-pass frequency characteristic with a cutoff and anti-aliasing filter before quantization. The amplifier also has non-linear distortion. The output is the signed 16 bit number (-2^15 to 2^15-1) where 2^15 corresponds to 3 volt amplitude on the quantizer.  This output is in an int16 (short integer) format and may need to be converted to a double format before applying other Matlab computations.

ampsqb.p - This function simulates digitizing a signal with nonlinear (amplifier) characteristics, a 16 bit quantization over a ± 3 volt range, and a sampling rate of 16kHz. The input should be a signal sampled at a rate higher than or equal to 16kHz for the simulation to be effective (if a lower rate is used, it will be upsampled, which is not realistic for an analog to digital converter).  The function is called by:

>>    sigout = ampsqb(sigin, fso)

where SIGIN represents the analog signal (not quantized) and created with a sampling rate of FSO Hz. The signal is passed through an anti-aliasing filter and sampled at 16kHz (resamples).  The amplifier has a high-pass frequency characteristic with a cutoff and anti-aliasing filter before quantization. The amplifier also has non-linear distortion. The output is the signed 16 bit number (-2^15 to 2^15-1) where 2^15 corresponds to 3 volt amplitude on the quantizer.  This output is in an int16 (short integer) format and may need to be converted to a double format before applying other Matlab computations.

modroom.p - This function inputs a source signal, its source position, and up to two receivers (microphones) at given positions.  The function syntax is given by:

>>  sigout = modroom(sigin,fs,sigpos, micpos1, micpos2)

where SIGIN is a vector containing the source signal sampled at rate FS and at [x, y, z] position denoted by vector SIGPOS in meters.  The microphone position is given by a similar vector for MICPOS1.  The output SIGOUT is the signal received at MICPOS1 with other noises and reverberations in the room, as well as quantization, distortion, and sampling to 16 kHz form the amp and quantizer pfile AMPSQA.  The augment MICPOS2 is optional.  If MICPOS2 is given, then SIGOUT will be a 2 column vector, resulting from signals received at both microphones. The simulated room recording includes additive noise signals from sources such as vents and buzzing from light and fans. Therefore if SIGIN is all zeros the output SIGOUT will consist of room noise only.  The room is rectangular with reflection coefficients relatively high for the 4 walls.  The reflection from the ceiling and floor is negligible.  The coordinate [0, 0, 0] is the center of the room. A cubic region of ±1.5 meters around the room center is a valid place to place sources and microphones. The actual walls and ceiling are beyond these points.

cclip.m  This function center clips a signal according to several option described in the function help.

featureex.m  This script presents an example of a Monte Carlo simulation for classifying segments containing either a simple tone with noise or noise only.  Only one feature is computed from the signal characterizations, the rest are left as an exercise.  This script requires the function tone.m to run.

mindistex.m  This script simulates a minimum distance classification problem where training data is generated for a 2 feature classifier and minimum distance classifiers are estimated; one without scaling and one with scaling by the covariance matrix.  Distances are referred to Mahanolobis distances when scaled by the covariance matrix.  Test data are independently generated to estimate the classification error of the classifiers derived from the training data.  The parameters are set to give an approximate classification error of 4% without scaling and 1% when scaling with the Mahanolobis distance.

ldaex2.m  This script presents an example of a Monte Carlo simulation for computing a linear discriminate for 2 classes based on mean and standard deviation statistics for the feature vectors.  The user must type these statistics in for each element in the feature vector.  The script plots out the feature space with projection and decision boundaries for 2-D feature vectors.  The script requires roc.m to compute the classification error for each Monte Carlo run.  In the end it presents a plot of the percent classification error for each run.

ldaexa.m  This script presents an example of a Monte Carlo simulation for computing a linear discriminate for 2 classes based on mean and standard deviation statistics for the feature vectors.  For each class, the user types in a single mean and variance for all features in the vector.  This is convenient for testing the effects of increasing the number of features.  The script plots out the feature space with projection and decision boundaries for 2-D feature vectors.  The script requires roc.m to compute the classification error for each Monte Carlo run.  In the end it presents a plot of the percent classification error for each run.

peakfind.m  This function inputs a sampled waveform and corresponding x-axis and returns the locations and amplitudes of all local maxima (it always include endpoints).

lpcex.m  This script presents an example of computing LPC coefficients for a sample of voiced speech (needs data file  aaa3.wav to run).  In addition, the root of the LPC are computed and used as an estimate of the formant frequencies.

lpctasks.m  This script presents an example of computing LPC coefficients for a sample of voiced speech (needs data file  aaa1.wav to run).  In addition, the estimated formants are used to estimate the length of the vocal fold track, reprocess/filter the signal with LPC filter and its inverse, and manipulate pole placements for the reconstructed sound.

ceppitchest.m  This function applies Noll's cepstrum pitch detection algorithm to an input signal vector.  Special processing parameters can be set through an input data structure.  If the input signal is longer than the set frame length, the function estimates the pitch for each frame as it is slid over the entire signal. The function returns vector containing the resulting pitch contour vector, with a corresponding detection statistic vector (i.e. numbers proportional to how likely the segment contained a signal with an actual pitch), and time value vector corresponding to each estimate.  This function calls another functions cclip.m; however, it is included at the end of this mfile and does not need to be downloaded separately.

paravoiced.m  This function estimates useful parameters for analyzing voiced signal segments.  This input is a data structure with the signal segment and sampling rate.  The output has various parameters including pitch and formant frequencies.  Example voiced speech segments can be downloaded from aah.wav, ooo.wav, and eee.wav.  This program requires ceppitchest.m to run.