Signal Processing Examples
The SignalProcessing package features tools for frequency domain analysis, windowing, signal generation and analysis, and more.
Getting Started
While any command in the package can be referred to using the long form, for example, SignalProcessing:-Spectrogram, it is often easier to load the package and then use the short form command names.
restart;
with(SignalProcessing):
Spectrogram, Power Spectrum and Signal Plots
Imaging a Speech Spectrogram
This is a spectrogram, power spectrum and signal plot of a male voice saying MapleSim, recorded at 11 kHz.
maplesim := FileTools:-JoinPath([ kernelopts(datadir), "audio", "maplesim.wav" ]):
Spectrogram(maplesim, colorscheme = ["zgradient", [white, LightSteelBlue, black], markers = [0, .5, 1]], fftsize = 256, includesignal, includepowerspectrum);
Spectrogram of a Violin Note
This is a spectrogram of a violin note played with vibrato, together with the power spectrum and signal plot. Note that the frequency oscillations are at 7 kHz and above, which is typical of vibrato.
filename := FileTools:-JoinPath([ kernelopts(datadir), "audio", "ViolinThreePosVibrato.wav" ]):
vibrato := AudioTools:-Read(filename):
Spectrogram(vibrato, colorscheme = ["zgradient", [white, LightBlue, red, black], markers = [0, .5, .75, 1]], channel = 1, fftsize = 2^10, includepowerspectrum, includesignal);
Filtering Signals
filename2 := FileTools:-JoinPath([ kernelopts(datadir), "audio", "MapleSimMono11025.wav" ]):
originalSpeech := AudioTools:-Read(filename2):
Plot waveform and power spectrum:
commonPlotOpts2 := labeldirections = [horizontal, vertical], axis = [gridlines = [color = "SteelBlue"]]:
samplingRate := attributes(originalSpeech)[1];
samplingRate≔11025
duration := evalf(AudioTools:-Duration(originalSpeech));
duration≔4.504671202
p1 := plots:-listplot([seq([i/samplingRate, originalSpeech[i]], i = 1 .. numelems(originalSpeech))], thickness = 0, gridlines, axes = boxed, title = "Original Speech", labels = ["Time (s)", Waveform], commonPlotOpts2):
plots:-display(p1);
fq := FFT(originalSpeech[1 .. 2^15]):
psq := PowerSpectrum(fq):
ps1 := plots:-pointplot([seq([i*samplingRate/2^15, psq[i]], i = 1 .. (1/2)*2^15)], thickness = 0, color = black, gridlines, connect = true, title = "Power Spectrum of Original Speech", labels = ["Frequency (Hz)", Power], view = [100 .. 2000, 0 .. 1.6], axis[1] = [mode = log], commonPlotOpts2):
plots:-display(ps1);
Apply IIR Butterworth or Chebyshev Filter:
fc := 800:
taps := GenerateButterworthTaps(9, fc/samplingRate, 'filtertype' = 'lowpass', 'normalize' = true):
filteredSpeech := InfiniteImpulseResponseFilter(originalSpeech, taps):
View before and after power spectrum and waveform:
FFTfilteredSpeech := FFT(filteredSpeech[1 .. 2^15]):
PSfilteredSpeech := PowerSpectrum(FFTfilteredSpeech):
ps2 := plots:-pointplot([seq([i*samplingRate/2^15, PSfilteredSpeech[i]], i = 1 .. (1/2)*2^15)], thickness = 0, color = black, gridlines, connect = true, title = "Power Spectrum of Filtered Speech", labels = ["Frequency (Hz)", Power], view = [100 .. 2000, 0 .. 1.6], axis[1] = [mode = log], commonPlotOpts2):
plots:-display(Array(<<ps1, ps2>>));
p2 := plots:-listplot([seq([i/samplingRate, filteredSpeech[i]], i = 1 .. numelems(originalSpeech))], thickness = 0, gridlines, axes = boxed, color = black, title = "Filtered Speech", labels = ["Time (s)", Waveform], commonPlotOpts2):
plots:-display(Array(<<p1, p2>>), view = [0 .. duration, -1 .. 1]);
Apply FIR filter
flow := 200:
fhigh := 700:
taps := GenerateFiniteImpulseResponseFilterTaps(50, [flow/samplingRate, fhigh/samplingRate], filtertype = bandpass):
filteredSpeech := FiniteImpulseResponseFilter(originalSpeech, taps):
Windowing Functions
Windowing functions can also be applied to each window of data used to generate an FFT.
freq := 5:
N := 2^8:
samplingRate := 50:
data := Vector(N, proc (i) options operator, arrow; sin(2*evalf(Pi)*freq*i/samplingRate)+(1/1000000000000)*rand() end proc, datatype = float[8]):
Apply a Blackman Nuttall Window:
dataWindow := BlackmanNuttallWindow(data):
Spectrogram(dataWindow, samplerate = samplingRate, includepowerspectrum, includesignal, fftsize = 128, colorscheme = ["zgradient", [SteelBlue, PaleGreen, red, Black], markers = [0, .5, .8, 1]]);
Visualizations using Explore
It is possible to create interactive applications that use commands from the SignalProcessing package using the Explore command. To see another example of using Explore to create For more examples, see the Filtering Frequency Domain Noise application. The following example dynamically illustrates the effects of a filter parameter value on the resulting modified signal.
Generate a known signal:
A := SignalProcessing:-GenerateJaehne( 512, 4095 ):
PA := plots:-listplot( A, 'title' = "Original Signal" ):
Next, construct a reusable container Array C in order to improve performance of the exploration process by leveraging the in-place computational semantics of the InfiniteImpulseResponseFilter command.
C := Array( 1..512, datatype=float[8] ):
F := proc( frequency, filter ) uses SignalProcessing; local PB,taps; taps:=GenerateButterworthTaps( 3, frequency, 'filtertype' = filter ); InfiniteImpulseResponseFilter( A, taps, 'container' = C ); PB:=plots:-listplot( C,'title'="Filtered Signal" ); plots:-display( <<PA;PB>>, 'size'=[700,400] ); end proc:
Explore( F( frequency, filter ), parameters=[frequency = 0.05..0.49, [filter = ["highpass","lowpass"]]], placement=left, animate, loop, title="Infinite Impulse Response Filter");
frequency
filter
highpasslowpass
See Also
applications/SignalGeneration, applications/SunspotPeriodicity, SignalProcessing, SignalProcessing[FFT], SignalProcessing[FiniteImpulseResponseFilter], SignalProcessing[GenerateButterworthTaps], SignalProcessing[GenerateFiniteImpulseResponseFilterTaps], SignalProcessing[InfiniteImpulseResponseFilter], SignalProcessing[PowerSpectrum]
Download Help Document
