Signal and Image Processing - New Features in Maple 2019 - Maplesoft

What's New in Maple 2019

Signal and Image Processing



Frequency Spectrum of Irregularly Sampled Data 

The new LSPeriodogram command uses the Lomb-Scargle technique to generate a periodogram of data sampled at irregular time intervals. This is in contrast to the standard Fourier approach (such as that used in Periodogram) which requires regularly sampled data.
 

> with(SignalProcessing); -1
 

Generate an irregularly spaced vector of times 

> t := sort(LinearAlgebra:-RandomVector(`^`(2, 10), generator = 0 .. 36.0, datatype = float[8]))
 

Typesetting:-msub(Typesetting:-mi(
 

Generate a signal using frequencies of 1 Hz and 8 Hz 

> f1 := 1.0; -1; f2 := 8.0; -1
 

> s := Vector(`^`(2, 10), proc (i) options operator, arrow; `+`(sin(`+`(`*`(2, `*`(f1, `*`(Pi, `*`(t[i])))))), `*`(1.5, `*`(sin(`+`(`*`(2, `*`(f2, `*`(Pi, `*`(t[i]))))))))) end proc, datatype = float[8]...
 

Typesetting:-msub(Typesetting:-mi(
 

Generate a periodogram assuming that the signal is sampled at each point in the irregularly spaced time vector 

> LSPeriodogram(t, s, frequencyscale =
 

Plot_2d
 

Several options allow you to specify the frequencies over which the periodogram is generated, the normalization and more. Additionally, you can also generate the numerical data used to create this plot with the LSSpectrum command. 

> LSSpectrum(t, s)
 

 

 

FindPeakPoints 

The new FindPeakPoints command finds the peaks and valleys of 1D data sets. This versatile command offers several options that let you filter out peaks or valleys that are too close, what defines a peak or valley, and more. 

Here, we find the fundamental frequency and harmonics of a violin note by locating the peak points of its periodogram. This speaker component is needed to play audio: Embedded component 

First import, play the data, and view the power spectrum. 

> violinNote := AudioTools:-Read(FileTools:-JoinPath([kernelopts(datadir),
violinNote := AudioTools:-Read(FileTools:-JoinPath([kernelopts(datadir),
violinNote := AudioTools:-Read(FileTools:-JoinPath([kernelopts(datadir),
 

Typesetting:-mprintslash([violinNote := MATRIX([[
 

> violinPeriodogram := Periodogram(violinNote, powerscale = plots:-display(violinPeriodogram, size = [800, 300], view = [0 .. 10000, default])
 

Plot_2d
 

Now isolate the peak points on the periodogram, while filtering out any spurious peaks. 

> periodogramData := plottools:-getdata(violinPeriodogram)[3]; -1peakPoints := FindPeakPoints(periodogramData, minimumheight = -40, minimumprominence = 40, includeendpoints = false)
 

Typesetting:-msub(Typesetting:-mi(
 

Overlay the peak points over the periodogram 

> peakPlot := plot(peakPoints, style = point, symbol = solidcircle, symbolsize = 15); -1plots:-display(peakPlot, violinPeriodogram, view = [1 .. 10000, default], size = [800, 300])
 

Plot_2d
 

The fundamental frequency of the violin note is the first frequency in peakPoints (987.8 Hz). The other values are the harmonics (which are approximately integer multiples of the fundamental frequency) 

> `~`[`/`](peakPoints[2 .. (), 1], peakPoints[1, 1])
 

Typesetting:-msub(Typesetting:-mi(
 

Spectrogram 

You can now specify the overlap when generating a spectrogram (in prior releases, the overlap was fixed at 50%). 

Increasing overlap increases frequency resolution, but at the expense of lower time resolution. In other words, increasing overlap will reveal smaller time events, but will smear those events over time. Increasing overlap also increases computation time because there are more time slices to analyze. 

 

Here, we import an audio file and generate a spectrogram with overlaps of 10% and 90%. 

 

> aud := AudioTools:-Read(FileTools:-JoinPath([kernelopts(datadir),
 

Typesetting:-mprintslash([aud := MATRIX([[
 

> Spectrogram(aud, overlap = .1, fftsize = 256, size = [800, 300])
 

Plot_2d
 

> Spectrogram(aud, overlap = .90, fftsize = 256, size = [800, 300])
 

Plot_2d
 

 

Real and Complex Cepstrum 

Maple 2019 features new commands for calculation the real cepstrum, complex cepstrum and inverse complex cepstrum of a signal. These have many applications in seismic analysis, the characterisation of explosions, speech analysis, homomorphic filtering, and more. 

In this example, the location of an echo in an audio file is found using the RealCepstrum command; this information is used to remove the echo with an IIR filter.  

> echoAudio := AudioTools:-Read(FileTools:-JoinPath([kernelopts(datadir), Fs := attributes(echoAudio)[1]
 

 

Typesetting:-mprintslash([echoAudio := MATRIX([[
Typesetting:-mprintslash([Fs := 22050], [22050])
 

Play the audio - note the echo 

> DocumentTools:-SetProperty(AudioTools:-Play(echoAudio,
 

 

Vector of times for each point in the audio signal 

> t := Vector(numelems(echoAudio), proc (i) options operator, arrow; `/`(`*`(1.0, `*`(i)), `*`(Fs)) end proc); -1
 

 

Use RealCepstrum to identify start of echo 

> c := RealCepstrum(echoAudio); -1plot(t, c, view = [default, -.3 .. .3], labels = [
 

Plot_2d
 

 

 

The cepstrum plot shows a strong peak at a quefrency of just over 0.4 s - that's probably where the echo starts. Let's find the precise location of the peak. 

> threshold := map(proc (x) options operator, arrow; piecewise(`<`(x, .1), 0, 1) end proc, c); -1ind := ArrayTools:-SearchArray(threshold, 4)
 

Typesetting:-mfenced(Typesetting:-mrow(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn(
 

The peak occurs at an index of 9042. Use this information to attenuate the echo with an IIR filter 

> filterCoeffs := Array([1, `$`(0, 9040), .5]); -1; cleanAudio := Filter(echoAudio, filterCoeffs, Array([1])); -1
 

 

Playing the audio reveals that the echo has been attenuated 

> AudioTools:-Play(cleanAudio,
 


Visualize waveform of audio before and after filtering
 

> pltEchoAudio := plot(t, echoAudio, color = pltCleanAudio := plot(t, cleanAudio, color = plots:-display(pltEchoAudio, pltCleanAudio)
 

Plot_2d
 

FFTShift 

FFTShift swaps data in a matrix or a vector into a different position. For a matrix swapped about both dimensions, the quadrants are moved like so. 

Image 

The Fourier transform of an image places lower frequency data near all four corners, with higher frequency data near the center. In this instance, FFTShift is typically applied to move the lowest frequencies to the center and the highest frequencies to the corners. 

This results in a more meaningful visualization of the power spectrum where the lowest frequency data is contiguous, and simplifies the manipulation of frequency data 

 

First, we import an image. 

> with(ImageTools); -1img := Read(FileTools:-JoinPath([kernelopts(datadir), Embed(img)
 

Image 

Now, we convert from the image domain to the spatial frequency domain, move the zero-frequency data to the center, and visualize the power spectrum (the square-root scaling simply rescales the frequencies for a more meaningful visualization). 

> img_fourier := FFTShift(FFT(img)); -1Embed(Create(`~`[sqrt](`~`[abs](img_fourier))))
 

Image 

Multiply the spatial frequencies by a Gaussian mask (this is similar to a low pass filter that attenuates higher frequency data). 

> nRows, nCols := LinearAlgebra:-Dimension(img); -1img_fourier_blur := `~`[`*`](img_fourier, ` $`, Matrix(nRows, nCols, proc (i, j) options operator, arrow; evalf(exp(`/`(`*`(-1, `*`(`+`(`*`(`^`(`+`(i, `-`(`*`(`/`(1, 2), `*`(nRows)))), 2)), `*`(`^`(`+...Embed(Create(`~`[sqrt](`~`[abs](img_fourier_blur))))
 

Image 

Move the zero frequency data back to their original positions, and then invert the spatial frequency data back to the image domain. 

> img_low_pass := `~`[Re](InverseFFT(FFTShift(img_fourier_blur))); -1Embed(FitIntensity(Create(img_low_pass)))
 

Image 

 

Edge Detect 

The new EdgeDetect command applies an edge detection convolution mask to an image. The command supports several masks including Sobel, Robert, Prewitt 3x3 and Prewitt 4x4. 

> edge1 := EdgeDetect(img, method = Sobel); -1; Embed(FitIntensity(edge1))
 

Image 

> edge2 := EdgeDetect(img, method = Prewitt4); -1; Embed(FitIntensity(edge2))
 

Image 

 

ImagePeriodogram 

ImagePeriodogram calculates the magnitude of the Fourier transform of an image; this data can then be embedded as an image to visualize the spatial frequencies of the image. By default the lowest frequencies are shifted to the center (this can, however, be modified). Moreover, you can also scale the data to better visualize spatial frequencies that would otherwise not be apparent). 

> img := Read(cat(kernelopts(datadir),
 

Image 

> Embed(ImagePeriodogram(img))
 

Image 

More Updates 

Detrend 

A new option for the SignalPlot, Periodogram and Spectrogram lets you remove a linear trend from the data, prior to the primary analysis.