Signal and Image Processing

All Products Maple MapleSim

Home : Support : Online Help : System : Information : Updates : Maple 2021 : Signal and Image Processing

Signal and Image Processing

The SignalProcessing and ImageTools packages have been expanded with new and updated commands along with enhanced tools in the Context Panel.

>	with( ImageTools ):

>	with( SignalProcessing ):

SignalProcessing

BandPower, MeanFrequency, and SpectralEntropy

The new BandPower, MeanFrequency, and SpectralEntropy commands are used to find the power spectral density of a signal, and compute, respectively, the band power, mean frequency, and spectral entropy, either for the entire signal, or a specific frequency band. For example:

>	sample_rate := 5000;

$sample_rate ≔ 5000$

(1.1.1)

>	Times := Vector( [ seq ]( 0 .. 1, 1.0 / sample_rate ), datatype = float[8] ):

>	g := t -> cos( 1000 * Pi * t ) + 3 * cos( 2000 * Pi * t ) + 2 * cos( 3000 * Pi * t ); Signal := Vector( g~(Times), datatype = float[8] ):

$g ≔ t \mapsto \cos (1000 \cdot π \cdot t) + 3 \cdot \cos (2000 \cdot π \cdot t) + 2 \cdot \cos (3000 \cdot π \cdot t)$

(1.1.2)

>	Periodogram( Signal, samplerate = sample_rate, size = [800,400] );

>	frequency_range := 250 * Unit( Hz ) .. 1.25 * Unit( kHz );

$frequency_range ≔ 250 ⟦Hz⟧ .. 1.25 ⟦kHz⟧$

(1.1.3)

>	BandPower( Signal, sample_rate, frequency_range ); MeanFrequency( Signal, sample_rate, frequency_range ); SpectralEntropy( Signal, sample_rate, frequency_range );

$5.00688167653967309$

$949.412589381238149$

$1.34562172301835759$

(1.1.4)

The commands can also return the power spectral density and frequencies:

>	Power_Spectral_Density, Frequencies := BandPower( Signal, sample_rate, output = [ psd, frequencies ] );

PowerSpectrum

The PowerSpectrum command now accepts signals (in addition to FFTs) and other new options. For instance:

>	n := 8192;

$n ≔ 8192$

(1.2.1)

>	fs := ( n - 1 ) / 2.0 / Pi;

$fs ≔ 1303.638139$

(1.2.2)

>	T := Vector( n, k -> (k-1) / fs, datatype = float[8] ):

>	X := Vector( n, k -> 3 * sin( 50 * T[k] ) + 5 * I * cos( 150 * T[k] ) + 10, datatype = complex[8] ):

>	R := PowerSpectrum( X, samplerate = fs, variety = signal, window = Hamming, powerscale = dB/Hz, output = record ):

>	'power' = R[ power ];

>	R[ periodogram ];

Welch

The Welch command estimates the power spectrum of a signal, while attenuating the effect of noise at the expense of frequency resolution. To this end, the signal is divided into overlapping segments of equal (or nearly equal) length, and for each segment, a window is applied and the power spectrum found. The overall power spectrum is determined by averaging all the segment power spectra.

For an example, consider the following signal:

>	num_points := 8192;

$num_points ≔ 8192$

(1.3.1)

>	sample_rate := ( n - 1 ) / 2.0 / Pi;

$sample_rate ≔ 1303.638139$

(1.3.2)

>	Times := Vector( num_points, k -> (k-1) / sample_rate, datatype = float[8] ):

>	g := t -> cos( 10 * Pi * t ) + 3 * cos( 20 * Pi * t ) + 2 * cos( 60 * Pi * t ); Original_Signal := Vector[column]( g~( Times ), datatype = float[8] ):

$g ≔ t \mapsto \cos (10 \cdot π \cdot t) + 3 \cdot \cos (20 \cdot π \cdot t) + 2 \cdot \cos (60 \cdot π \cdot t)$

(1.3.3)

Now, add noise:

>	use Statistics in Noise := Vector[column]( Sample( RandomVariable( Normal( 0, 2 ) ), num_points ) ); end use:

>	Noisy_Signal := Original_Signal + Noise:

Now, let's compare plots of the signals and power spectra:

>	plots:-display( Array( [ dataplot( Times, Original_Signal, style = line, view = [0..2Pi,-40..40], title = "Original Signal" ), dataplot( Times, Noisy_Signal, style = line, view = [0..2Pi,-40..40], title = "Noisy Signal" ) ] ) );

plots:-display( Array( [
Welch( Original_Signal, samplerate = sample_rate, segmentsize = num_points, output = periodogram, periodogramoptions = [ title = "Periodogram of Original Signal" ] ),
Welch( Noisy_Signal, samplerate = sample_rate, segmentsize = num_points, output = periodogram, periodogramoptions = [ title = "Periodogram of Noisy Signal" ] )
] ) );

Now, apply Welch's method with an appropriate choice of parameter values to return a record with the Vectors for the power spectrum and frequencies, along with the periodogram:

>	R := Welch( Noisy_Signal, overlapsize = 512, segmentsize = 1024, window = "Hamming", samplerate = sample_rate, temperendpoints = true, datapowerscale = 1/Hz, plotpowerscale = dB/rad/Hz, output = record ):

>	'power' = R[power];

>	'frequencies' = R[frequencies];

>	R[periodogram];

MUSIC

The MUSIC command performs the Multiple Signal Classifier (MUSIC) Method on a signal, which estimates frequencies present in a noisy signal by computing the eigenvalues of the autocorrelation matrix and separating them into signal-subspace eigenvalues and noise-subspace eigenvalues. Consider, for an example, the following real-valued signal:

>	# number of points in the signal m := 2^8:

>	# sample rate fs := 100.0:

>	# times T := Vector( m, k -> 2 * Pi * (k-1) / fs, 'datatype' = 'float[8]' ):

>	# pure signal X := Vector( m, k -> 5.00 * sin( 10.25 * T[k] ) + 3.00 * sin( 10.40 * T[k] ) - 7.00 * sin( 20.35 * T[k] ), 'datatype' = 'float[8]' ):

To this, we add some noise:

>	# noise use Statistics in N := Vector[column]( Sample( RandomVariable( Normal( 0, 1 ) ), m ) ): end use:

>	# noisy signal Y := X + N:

Let's compare plots of the signals and power spectra for both the original and noisy versions:

>	DocumentTools:-Tabulate( Array( [ SignalPlot( X, 'samplerate' = fs, 'color' = 'blue', 'title' = "Pure Signal Plot", 'font' = [ 'Verdana', 15 ] ), SignalPlot( Y, 'samplerate' = fs, 'color' = 'blue', 'title' = "Noisy Signal Plot", 'font' = [ 'Verdana', 15 ] ) ] ) ):

>	DocumentTools:-Tabulate( Array( [ Periodogram( X, 'samplerate' = fs, 'color' = 'blue', 'title' = "Pure Signal Periodogram", 'font' = [ 'Verdana', 15 ] ), Periodogram( Y, 'samplerate' = fs, 'color' = 'blue', 'title' = "Noisy Signal Periodogram", 'font' = [ 'Verdana', 15 ] ) ] ) ):

The power spectra for both the clean and noisy signals fails to discriminate the two smaller frequencies. With the MUSIC command, however, we can detect all the frequencies:

>	MUSIC( Y, 'samplerate' = fs, 'dimension' = 6, 'output' = 'plot' );

The peaks can also be returned, with the dominant peaks corresponding to the frequencies of the sinusoids in the clean signal:

>	MUSIC( Y, 'samplerate' = fs, 'dimension' = 6, 'output' = 'peaks' );

ShortTimeFourierTransform

The new ShortTimeFourierTransform command takes a signal, and, similar to the Welch command, divides it into overlapping segments of equal (or nearly equal) length, and for each segment, a window is applied and the Fourier Transform computed. For example, here we find the Short-Time Fourier Transform of a violin recording along with the spectrogram (which plots the corresponding Short-Time Power Spectrum versus time):

>	with( AudioTools ):

>	file := cat( kernelopts( datadir ), kernelopts( dirsep ), "audio", kernelopts( dirsep ), "ViolinThreePosVibrato.wav" );

$file ≔ C:\Program Files\Maple Main\data\audio\ViolinThreePosVibrato.wav$

(1.5.1)

>	violin := ToMono( Read( file, samples = 1000 .. 10000 ) );

$violin ≔ [\begin{array}{c} Sample Rate & 44100 \\ File Format & PCM \\ File Bit Depth & 16 \\ Channels & 1 \\ Samples/Channel & 9001 \\ Duration & 0.20410 s \end{array}]$

(1.5.2)

>	Periodogram( violin, size = [800,400] );

>	overlap_size := 128;

$overlap_size ≔ 128$

(1.5.3)

>	segment_size := 256;

$segment_size ≔ 256$

(1.5.4)

>	stft_matrix, spectrogram_plot := ShortTimeFourierTransform( violin, overlapsize = overlap_size, segmentsize = segment_size, frequencyunit = kHz, powerscale = dB/Hz, output = [ stft, spectrogram ] ):

>	stft_matrix;

>	spectrogram_plot;

ShortTimeBandPower, ShortTimeMeanFrequency, and ShortTimeSpectralEntropy

The new commands ShortTimeBandPower, ShortTimeMeanFrequency, and ShortTimeSpectralEntropy apply the ShortTimeFourierTransform command, and compute the respective statistics for each short-time interval. Continuing the example above for the violin, we use the ShortTimeBandPower command to return everything, including plots, in a record:

>	R := ShortTimeBandPower( violin, overlapsize = overlap_size, segmentsize = segment_size, frequencyunit = kHz, output = record ):

>	R[bandpowerplot];

>	R[meanfrequencyplot];

>	R[entropyplot];

FilterFrequencyResponse

The new FilterFrequencyResponse command takes one (for an FIR filter) or two (for an IIR filter) containers for taps, representing the coefficients of the transfer function, and returns the frequency response. For example:

>	T := GenerateChebyshev1Taps( 5, 0.4, 29, filtertype = lowpass );

$T ≔ [\begin{array}{c} 0.0145256222951220 & 0.0726281114756102 & 0.145256222951220 & 0.145256222951220 & 0.0726281114756102 & 0.0145256222951220 & 1. & 1.67764328574172 & 1.12542851668154 & −0.885659875792078 & −1.51371985423053 & −0.938872158956746 \end{array}]$

(1.7.1)

>	A, B := ListTools:-Slice( T, 2 );

$A, B ≔ [\begin{array}{c} 0.0145256222951220 & 0.0726281114756102 & 0.145256222951220 & 0.145256222951220 & 0.0726281114756102 & 0.0145256222951220 \end{array}], [\begin{array}{c} 1. & 1.67764328574172 & 1.12542851668154 & −0.885659875792078 & −1.51371985423053 & −0.938872158956746 \end{array}]$

(1.7.2)

>	R := FilterFrequencyResponse( A, B, output = record ):

>	'response' = R[ response ];

>	R[ magnitudeplot ];

>	R[ phaseplot ];

EquivalentNoiseBandwidth

The EquivalentNoiseBandwidth command computes the equivalent-noise bandwidth of a window. For instance:

>	EquivalentNoiseBandwidth( 10, ["Exponential",0.25], 1.5 );

$0.173009806903817592$

(1.8.1)

Hampel

The new Hampel command is useful when removing outliers from data. For an example, we will create a signal, add outliers, and then apply the filter:

>	numpoints := 51;

$numpoints ≔ 51$

(1.9.1)

>	SignalOriginal := Vector( numpoints, i -> 5 + cos( 4 * Pi * (i-1) / (numpoints-1) ), datatype = float[8] ):

>	SignalModified := copy( SignalOriginal ): SignalModified[3] := SignalOriginal[3] + 4.0: SignalModified[floor(numpoints/2)] := SignalOriginal[floor(numpoints/2)] + 2.5: SignalModified[-2] := SignalModified[-2] - 3.0: 'SignalModified' = SignalModified:

>	SignalFiltered := Hampel( SignalModified, 3, 2.0, output = hampelsignal );

dataplot(
    [ SignalOriginal, SignalModified, SignalFiltered ],
    view = [ 1 .. numpoints, 0 .. 10 ],
    labels = [time,amplitude],
    title = "Signals",
    font = [Verdana,15],
    legend = ["original signal","modified signal","filtered modified signal"],
    legendstyle = [font=[Verdana,15]],
    symbolsize = 5,
    color = ["green","red","blue"],
    thickness = [6,3,3]
);

IntegrateData

The new IntegrateData command is used to estimate the area beneath a 1-D signal. For the example below, we start with a symbolic expression to create a signal, so that we can compare the symbolic integral with the numeric approximation:

>	signal := t -> 5 + sin( t ) + 0.1 * sin( 3 * t ) + 0.01 * sin( 5 * t );

$signal ≔ t \mapsto 5 + \sin (t) + 0.1 \cdot \sin (3 \cdot t) + 0.01 \cdot \sin (5 \cdot t)$

(1.10.1)

>	t1, t2 := 0.0, 10.0;

$t1, t2 ≔ 0., 10.0$

(1.10.2)

>	points := 100;

$points ≔ 100$

(1.10.3)

>	tValues := Vector( [ seq( t1 .. t2, numelems = points ) ], datatype = float[8] ):

>	yValues := signal~( tValues ):

>	area_symbolic := int( signal, t1 .. t2 );

$area_symbolic ≔ 51.86733322$

(1.10.4)

>	area_numeric := IntegrateData( tValues, yValues, method = simpson );

$area_numeric ≔ 51.8673373335523635$

(1.10.5)

FindPeakPoints

The FindPeakPoints command has been updated to include a new calling sequence, FindPeakPoints(X,Y,...), and two new options, maximumheight and sortdata (which, when false, is used to skip sorting if the independent data are known to be sorted).

For example, consider the following signal:

>	f := unapply( sin(t) + 3 * cos(3t), t ); a, b := 0, 6 Pi; plot( f, a .. b, 'color' = 'blue' );

$f ≔ t \mapsto \sin (t) + 3 \cdot \cos (3 \cdot t)$

$a, b ≔ 0, 6 π$

We will create data from this, and add noise:

>	n := 10^4; T := Array( 1 .. n, i -> evalhf( (b-a) * (i-1) / (n-1) ), 'datatype' = 'float[8]' ); X := map['evalhf']( f, T ); use Statistics in N := Array( Sample( RandomVariable( Normal( 0, 0.1 ) ), n ) ); end use;

$n ≔ 10000$

The original signal has peaks of three different heights, and valleys of three different heights. Here, let's identify the smaller peaks and larger valleys in the noisy data:

>	FindPeakPoints( T, X + N, 'sortdata' = 'false', 'maximumheight' = 2.5, 'minimumheight' = -2.5, 'minimumprominence' = 0.75, 'output' = 'plot', 'plotincludepoints' = ['peaks','valleys'], 'font' = ['Verdana',15] );

ComplexToReal and RealToComplex

The new ComplexToReal and RealToComplex commands, respectively, combine containers with the real and imaginary parts into a single container with the complex values, and decompose a container with complex values into separate containers with the real and imaginary parts. For instance:

>	RealParts := LinearAlgebra:-RandomVector( 3, datatype = float[8] ); ImaginaryParts := LinearAlgebra:-RandomVector( 3, datatype = float[8] );

$RealParts ≔ [\begin{array}{c} 4. \\ 92. \\ −17. \end{array}]$

$ImaginaryParts ≔ [\begin{array}{c} 59. \\ −20. \\ −99. \end{array}]$

(1.12.1)

>	ComplexValues := RealToComplex( RealParts, ImaginaryParts );

$ComplexValues ≔ [\begin{array}{c} 4. + 59. I \\ 92. - 20. I \\ −17. - 99. I \end{array}]$

(1.12.2)

>	ComplexToReal( ComplexValues );

$[\begin{array}{c} 4. \\ 92. \\ −17. \end{array}], [\begin{array}{c} 59. \\ −20. \\ −99. \end{array}]$

(1.12.3)

The operations are optimized for hardware floats and large rtables, and existing containers can be passed for storage (using the container option for RealToComplex, and containers option for ComplexToReal).

RootMeanSquare

The RootMeanSquare command now supports multidimensional rtables and lists when computing the Root Mean Square (RMS). For instance:

>	M := Matrix( [ [ 1, 2 ], [ 3, 4 ] ], datatype = float[8] );

$M ≔ [\begin{array}{c} 1. & 2. \\ 3. & 4. \end{array}]$

(1.13.1)

>	RootMeanSquare( M );

$2.73861278752583059$

(1.13.2)

>	L := [ 10, 15, 20 ];

$L ≔ [10, 15, 20]$

(1.13.3)

>	RootMeanSquare( L );

$15.5456317551480261$

(1.13.4)

RootMeanSquareError and RelativeRootMeanSquareError

The new RootMeanSquareError and RelativeRootMeanSquareError commands are useful for quantifying errors. Specifically, the Root Mean Square Error (RMSE) of two containers X and Y is the RMS of the difference X-Y, and the Relative Root Mean Square Error (RRMSE) is the RMS of X-Y divided by the RMS of Y. For example:

>	RootMeanSquareError( < 1.1, 1.9, 3.1 >, [ 1, 2, 3 ] );

$0.100000000000000103$

(1.14.1)

>	P := < 0.00003, 0.00004 >; Q := < 0.00001, 0.00002 >;

$P ≔ [\begin{array}{c} 0.00003 \\ 0.00004 \end{array}]$

$Q ≔ [\begin{array}{c} 0.00001 \\ 0.00002 \end{array}]$

(1.14.2)

>	RootMeanSquareError( P, Q ); RelativeRootMeanSquareError( P, Q );

$0.0000199999999999999982$

$1.26491106406735154$

(1.14.3)

Mean

The Mean command in SignalProcessing now supports multidimensional containers and weights. For example:

>	A := Matrix( [ [ 5, 2 ], [ -1, 8 ] ], datatype = float[8] );

$A ≔ [\begin{array}{c} 5. & 2. \\ −1. & 8. \end{array}]$

(1.15.1)

>	Mean( A );

$3.50000000000000000$

(1.15.2)

>	W := Matrix( [ [ 1, 2 ], [ 3, 4 ] ], datatype = float[8] );

$W ≔ [\begin{array}{c} 1. & 2. \\ 3. & 4. \end{array}]$

(1.15.3)

>	Mean( A, W );

$3.79999999999999982$

(1.15.4)

Phase

The Phase command in SignalProcessing now includes an option for unwrapping the phases, so that there are no large jumps:

>	Z := Vector( 720, k -> k * exp( I * ( Pi / 180 * k )^2 ), datatype = complex[8] ):

>	P := Phase( Z ):

>	dataplot( P, style = line, title = "Wrapped Phase" );

>	Q := Phase( Z, unwrap );

>	dataplot( Q, style = line, title = "Unwrapped Phase" );

ImageTools

The SampleImage command in the ImageTools package returns the requested image from a repository of sample images. For example:

>	# Current number of available images in the repository. num_images := SampleImage( 0 );

$num_images ≔ 5$

(2.1)

>	Embed( SampleImage( 1 ) );

>	Embed( SampleImage( 5 ) );

Download Help Document

Maple

Maple Add-Ons

MapleSim

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim

Maple

Powerful math software that is easy to use

Maple Add-Ons

MapleSim

Advanced System Level Modeling

MapleSim Add-Ons

Systems Engineering

Consulting Services

Maple T.A. and Möbius

Education

Industries

Automotive and Aerospace

Robotics

Machine Design & Industrial Automation

Other

Application Areas

Product Pricing

Purchasing

Institutional Student Licensing

Maplesoft Elite Maintenance (EMP)

Support

Product Training

Online Product Help

Webinars & Events

Publications

Content Hubs

Examples & Applications

Community

About Maplesoft

Media Center

User Community

Contact

Online Help

All Products Maple MapleSim