Signal and Image Processing - Maple Help

Home : Support : Online Help : What's New and Release Notes : 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;
 ${\mathrm{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}{↦}{\mathrm{cos}}\left({1000}{\cdot }{\mathrm{\pi }}{\cdot }{t}\right){+}{3}{\cdot }{\mathrm{cos}}\left({2000}{\cdot }{\mathrm{\pi }}{\cdot }{t}\right){+}{2}{\cdot }{\mathrm{cos}}\left({3000}{\cdot }{\mathrm{\pi }}{\cdot }{t}\right)$ (1.1.2)
 > Periodogram( Signal, samplerate = sample_rate, size = [800,400] );
 > frequency_range := 250 * Unit( Hz ) .. 1.25 * Unit( kHz );
 ${\mathrm{frequency_range}}{≔}{250}{}⟦{\mathrm{Hz}}⟧{..}{1.25}{}⟦{\mathrm{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 ] );
 ${\mathrm{Power_Spectral_Density}}{,}{\mathrm{Frequencies}}{≔}\begin{array}{c}\left[\begin{array}{c}{2.87942720498115}{×}{{10}}^{{-6}}\\ {2.87943448146849}{×}{{10}}^{{-6}}\\ {2.87945631139094}{×}{{10}}^{{-6}}\\ {2.87949269545107}{×}{{10}}^{{-6}}\\ {2.87954363522815}{×}{{10}}^{{-6}}\\ {2.87960913240025}{×}{{10}}^{{-6}}\\ {2.87968918936491}{×}{{10}}^{{-6}}\\ {2.87978380934996}{×}{{10}}^{{-6}}\\ {2.87989299596818}{×}{{10}}^{{-6}}\\ {2.88001675264265}{×}{{10}}^{{-6}}\\ {2.88015508507797}{×}{{10}}^{{-6}}\\ {2.88030799748615}{×}{{10}}^{{-6}}\\ {2.88047549605390}{×}{{10}}^{{-6}}\\ {2.88065758716316}{×}{{10}}^{{-6}}\\ {2.88085427745574}{×}{{10}}^{{-6}}\\ {⋮}\end{array}\right]\\ \hfill {\text{2501 element Vector[column]}}\end{array}{,}\begin{array}{c}\left[\begin{array}{c}{0.}\\ {0.999800039992002}\\ {1.99960007998400}\\ {2.99940011997600}\\ {3.99920015996801}\\ {4.99900019996001}\\ {5.99880023995201}\\ {6.99860027994401}\\ {7.99840031993601}\\ {8.99820035992802}\\ {9.99800039992002}\\ {10.9978004399120}\\ {11.9976004799040}\\ {12.9974005198960}\\ {13.9972005598880}\\ {⋮}\end{array}\right]\\ \hfill {\text{2501 element Vector[column]}}\end{array}$ (1.1.5)

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;
 ${\mathrm{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 ];
 ${\mathrm{power}}{=}{{\mathrm{_rtable}}}_{{36893489685096521348}}$ (1.2.3)
 > 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;
 ${\mathrm{num_points}}{≔}{8192}$ (1.3.1)
 > sample_rate := ( n - 1 ) / 2.0 / Pi;
 ${\mathrm{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}{↦}{\mathrm{cos}}\left({10}{\cdot }{\mathrm{\pi }}{\cdot }{t}\right){+}{3}{\cdot }{\mathrm{cos}}\left({20}{\cdot }{\mathrm{\pi }}{\cdot }{t}\right){+}{2}{\cdot }{\mathrm{cos}}\left({60}{\cdot }{\mathrm{\pi }}{\cdot }{t}\right)$ (1.3.3)

 > 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..2*Pi,-40..40], title = "Original Signal" ),     dataplot( Times, Noisy_Signal, style = line, view = [0..2*Pi,-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];
 ${\mathrm{power}}{=}{{\mathrm{_rtable}}}_{{36893489685096545684}}$ (1.3.4)
 > 'frequencies' = R[frequencies];
 ${\mathrm{frequencies}}{=}{{\mathrm{_rtable}}}_{{36893489685096545804}}$ (1.3.5)
 > 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' );
 ${{\mathrm{_rtable}}}_{{36893489685096513276}}$ (1.4.1)

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" );
 ${\mathrm{file}}{≔}{"C:\Program Files\Maple Main\data\audio\ViolinThreePosVibrato.wav"}$ (1.5.1)
 > violin := ToMono( Read( file, samples = 1000 .. 10000 ) );
 ${\mathrm{violin}}{≔}\left[\begin{array}{cc}{"Sample Rate"}& {44100}\\ {"File Format"}& {\mathrm{PCM}}\\ {"File Bit Depth"}& {16}\\ {"Channels"}& {1}\\ {"Samples/Channel"}& {9001}\\ {"Duration"}& {0.20410}{}{s}\end{array}\right]$ (1.5.2)
 > Periodogram( violin, size = [800,400] );
 > overlap_size := 128;
 ${\mathrm{overlap_size}}{≔}{128}$ (1.5.3)
 > segment_size := 256;
 ${\mathrm{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;
 ${{\mathrm{_rtable}}}_{{36893489685212346476}}$ (1.5.5)
 > 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 );
 $\left[\begin{array}{cccccccccccc}0.014525622295122047& 0.07262811147561024& 0.14525622295122048& 0.14525622295122048& 0.07262811147561024& 0.014525622295122047& 1.0& 1.6776432857417185& 1.1254285166815392& -0.8856598757920775& -1.5137198542305283& -0.9388721589567464\end{array}\right]$ (1.7.1)
 > A, B := ListTools:-Slice( T, 2 );
 ${A}{,}{B}{≔}\left[\begin{array}{cccccc}{0.0145256222951220}& {0.0726281114756102}& {0.145256222951220}& {0.145256222951220}& {0.0726281114756102}& {0.0145256222951220}\end{array}\right]{,}\left[\begin{array}{cccccc}{1.}& {1.67764328574172}& {1.12542851668154}& {-0.885659875792078}& {-1.51371985423053}& {-0.938872158956746}\end{array}\right]$ (1.7.2)
 > R := FilterFrequencyResponse( A, B, output = record ):
 > 'response' = R[ response ];
 ${\mathrm{response}}{=}{{\mathrm{_rtable}}}_{{36893489685173134268}}$ (1.7.3)
 > 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;
 ${\mathrm{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 );
 ${{\mathrm{_rtable}}}_{{36893489685185931372}}$ (1.9.2)
 > 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 );
 ${\mathrm{signal}}{≔}{t}{↦}{5}{+}{\mathrm{sin}}\left({t}\right){+}{0.1}{\cdot }{\mathrm{sin}}\left({3}{\cdot }{t}\right){+}{0.01}{\cdot }{\mathrm{sin}}\left({5}{\cdot }{t}\right)$ (1.10.1)
 > t1, t2 := 0.0, 10.0;
 ${\mathrm{t1}}{,}{\mathrm{t2}}{≔}{0.}{,}{10.0}$ (1.10.2)
 > points := 100;
 ${\mathrm{points}}{≔}{100}$ (1.10.3)
 > tValues := Vector( [ seq( t1 .. t2, numelems = points ) ], datatype = float[8] ):
 > yValues := signal~( tValues ):
 > area_symbolic := int( signal, t1 .. t2 );
 ${\mathrm{area_symbolic}}{≔}{51.86733322}$ (1.10.4)
 > area_numeric := IntegrateData( tValues, yValues, method = simpson );
 ${\mathrm{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(3*t), t ); a, b := 0, 6 * Pi; plot( f, a .. b, 'color' = 'blue' );
 ${f}{≔}{t}{↦}{\mathrm{sin}}\left({t}\right){+}{3}{\cdot }{\mathrm{cos}}\left({3}{\cdot }{t}\right)$
 ${a}{,}{b}{≔}{0}{,}{6}{}{\mathrm{\pi }}$

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}$
 ${{\mathrm{_rtable}}}_{{36893489685185597060}}$
 ${{\mathrm{_rtable}}}_{{36893489685185591284}}$
 ${{\mathrm{_rtable}}}_{{36893489685096560748}}$ (1.11.1)

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] );
 $\left[\begin{array}{c}4.0\\ 92.0\\ -17.0\end{array}\right]$
 $\left[\begin{array}{c}59.0\\ -20.0\\ -99.0\end{array}\right]$ (1.12.1)
 > ComplexValues := RealToComplex( RealParts, ImaginaryParts );
 $\left[\begin{array}{c}4.0+59.0{}I\\ 92.0-20.0{}I\\ -17.0-99.0{}I\end{array}\right]$ (1.12.2)
 > ComplexToReal( ComplexValues );
 $\left[\begin{array}{c}{4.}\\ {92.}\\ {-17.}\end{array}\right]{,}\left[\begin{array}{c}{59.}\\ {-20.}\\ {-99.}\end{array}\right]$ (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] );
 $\left[\begin{array}{cc}1.0& 2.0\\ 3.0& 4.0\end{array}\right]$ (1.13.1)
 > RootMeanSquare( M );
 ${2.73861278752583059}$ (1.13.2)
 > L := [ 10, 15, 20 ];
 ${L}{≔}\left[{10}{,}{15}{,}{20}\right]$ (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 >;
 $\left[\begin{array}{c}0.00003\\ 0.00004\end{array}\right]$
 $\left[\begin{array}{c}0.00001\\ 0.00002\end{array}\right]$ (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] );
 $\left[\begin{array}{cc}5.0& 2.0\\ -1.0& 8.0\end{array}\right]$ (1.15.1)
 > Mean( A );
 ${3.50000000000000000}$ (1.15.2)
 > W := Matrix( [ [ 1, 2 ], [ 3, 4 ] ], datatype = float[8] );
 $\left[\begin{array}{cc}1.0& 2.0\\ 3.0& 4.0\end{array}\right]$ (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 );
 ${{\mathrm{_rtable}}}_{{36893489685163984044}}$ (1.16.1)
 > 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 );
 ${\mathrm{num_images}}{≔}{5}$ (2.1)
 > Embed( SampleImage( 1 ) );

 > Embed( SampleImage( 5 ) );