Signal Processing - Maple Help

 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.

 > $\mathrm{with}\left(\mathrm{SignalProcessing}\right):$

Generate an irregularly spaced vector of times

 >
 ${{\mathrm{_rtable}}}_{{18446746419362836598}}$ (1.1)

Generate a signal using frequencies of 1 Hz and 8 Hz

 > $\mathrm{f1}≔1.0:\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{f2}≔8.0:$
 > $s≔\mathrm{Vector}\left({2}^{10},i→\mathrm{sin}\left(\mathrm{f1}\cdot 2\mathrm{\pi }{t}_{i}\right)+1.5\mathrm{sin}\left(\mathrm{f2}\cdot 2\mathrm{\pi }{t}_{i}\right),\mathrm{datatype}={\mathrm{float}}_{8}\right)$
 ${{\mathrm{_rtable}}}_{{18446746419350756518}}$ (1.2)

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

 >

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.

 > $\mathrm{LSSpectrum}\left(t,s\right)$
 $\begin{array}{c}\left[\begin{array}{c}{0.00556035750021131}\\ {0.0111207150004226}\\ {0.0166810725006339}\\ {0.0222414300008452}\\ {0.0278017875010566}\\ {0.0333621450012679}\\ {0.0389225025014792}\\ {0.0444828600016905}\\ {0.0500432175019018}\\ {0.0556035750021131}\\ {⋮}\end{array}\right]\\ \hfill {\text{2560 element Vector[column]}}\end{array}{,}\begin{array}{c}\left[\begin{array}{c}{0.00698512344159072}\\ {0.00776423599903224}\\ {0.00716845005457644}\\ {0.00610279006191028}\\ {0.00499610189423538}\\ {0.00538757255816807}\\ {0.00732037232391206}\\ {0.00768212344624377}\\ {0.00573470249368044}\\ {0.00301752080536450}\\ {⋮}\end{array}\right]\\ \hfill {\text{2560 element Vector[column]}}\end{array}$ (1.3)

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:

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

 >
 ${\mathrm{violinNote}}{≔}\left[\begin{array}{cc}{"Sample Rate"}& {44100}\\ {"Bit Depth"}& {16}\\ {"Channels"}& {1}\\ {"Points/Channel"}& {64724}\\ {"Duration"}& {1.47}{}{s}\end{array}\right]$ (2.1)
 > $\mathrm{violinPeriodogram}≔\mathrm{Periodogram}\left(\mathrm{violinNote},\mathrm{powerscale}="dB"\right):\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}$$\mathrm{plots}:-\mathrm{display}\left(\mathrm{violinPeriodogram},\mathrm{size}=\left[800,300\right],\mathrm{view}=\left[0..10000,\mathrm{default}\right]\right)$

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

 > $\mathrm{periodogramData}≔{\mathrm{plottools}:-\mathrm{getdata}\left(\mathrm{violinPeriodogram}\right)}_{3}:$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{peakPoints}≔\mathrm{FindPeakPoints}\left(\mathrm{periodogramData},\mathrm{minimumheight}=-40,\mathrm{minimumprominence}=40,\mathrm{includeendpoints}=\mathrm{false}\right)$
 ${{\mathrm{_rtable}}}_{{18446746419590517934}}$ (2.2)

Overlay the peak points over the periodogram

 >

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)

 > $\mathrm{peakPoints}\left[2..,1\right]/~\mathrm{peakPoints}\left[1,1\right]$
 ${{\mathrm{_rtable}}}_{{18446746419590198566}}$ (2.3)

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%.

 >
 ${\mathrm{aud}}{≔}\left[\begin{array}{cc}{"Sample Rate"}& {11025}\\ {"Bit Depth"}& {16}\\ {"Channels"}& {1}\\ {"Points/Channel"}& {8227}\\ {"Duration"}& {0.75}{}{s}\end{array}\right]$ (3.1)
 > $\mathrm{Spectrogram}\left(\mathrm{aud},\mathrm{overlap}=0.1,\mathrm{fftsize}=256,\mathrm{size}=\left[800,300\right]\right)$
 > $\mathrm{Spectrogram}\left(\mathrm{aud},\mathrm{overlap}=0.90,\mathrm{fftsize}=256,\mathrm{size}=\left[800,300\right]\right)$

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.

 > $\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{Fs}≔{\mathrm{attributes}\left(\mathrm{echoAudio}\right)}_{1}$
 ${\mathrm{echoAudio}}{≔}\left[\begin{array}{cc}{"Sample Rate"}& {22050}\\ {"Bit Depth"}& {16}\\ {"Channels"}& {1}\\ {"Points/Channel"}& {29232}\\ {"Duration"}& {1.33}{}{s}\end{array}\right]$
 ${\mathrm{Fs}}{≔}{22050}$ (4.1)

Play the audio - note the echo

 > $\mathrm{DocumentTools}:-\mathrm{SetProperty}\left("Speaker0",\mathrm{samplerate},\mathrm{Fs}\right);$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{AudioTools}:-\mathrm{Play}\left(\mathrm{echoAudio},"Speaker0"\right)$

Vector of times for each point in the audio signal

 > $t≔\mathrm{Vector}\left(\mathrm{numelems}\left(\mathrm{echoAudio}\right),i→\frac{1.0i}{\mathrm{Fs}}\right):$

Use RealCepstrum to identify start of echo

 > $c≔\mathrm{RealCepstrum}\left(\mathrm{echoAudio}\right):$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{plot}\left(t,c,\mathrm{view}=\left[\mathrm{default},-0.3..0.3\right],\mathrm{labels}=\left["Quefrency \left(s\right)",""\right],\mathrm{labeldirections}=\left[\mathrm{horizontal},\mathrm{vertical}\right]\right)$

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.

 > $\mathrm{threshold}≔\mathrm{map}\left(x→\mathrm{piecewise}\left(x<0.1,0,1\right),c\right):\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}$$\mathrm{ind}≔\mathrm{ArrayTools}:-\mathrm{SearchArray}\left(\mathrm{threshold},4\right)$
 $\left[\begin{array}{r}2\\ 7\\ 29\\ 9042\end{array}\right]$ (4.2)

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

 >

Playing the audio reveals that the echo has been attenuated

 > $\mathrm{AudioTools}:-\mathrm{Play}\left(\mathrm{cleanAudio},"Speaker0"\right)$

Visualize waveform of audio before and after filtering

 > $\mathrm{pltEchoAudio}≔\mathrm{plot}\left(t,\mathrm{echoAudio},\mathrm{color}="DarkOliveGreen",\mathrm{legend}="Audio with echo"\right):\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}$$\mathrm{pltCleanAudio}≔\mathrm{plot}\left(t,\mathrm{cleanAudio},\mathrm{color}="black",\mathrm{legend}="Audio with echo removed"\right):$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{plots}:-\mathrm{display}\left(\mathrm{pltEchoAudio},\mathrm{pltCleanAudio}\right)$

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.

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.

 > $\mathrm{with}\left(\mathrm{ImageTools}\right):\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}$$\mathrm{Embed}\left(\mathrm{img}\right)$

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).

 > $\mathrm{img_fourier}≔\mathrm{FFTShift}\left(\mathrm{FFT}\left(\mathrm{img}\right)\right):$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{Embed}\left(\mathrm{Create}\left(\mathrm{sqrt}~\left(\mathrm{abs}~\left(\mathrm{img_fourier}\right)\right)\right)\right)$

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

 > $\mathrm{nRows},\mathrm{nCols}≔\mathrm{LinearAlgebra}:-\mathrm{Dimension}\left(\mathrm{img}\right):\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{Embed}\left(\mathrm{Create}\left(\mathrm{sqrt}~\left(\mathrm{abs}~\left(\mathrm{img_fourier_blur}\right)\right)\right)\right)$

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

 > $\mathrm{img_low_pass}≔\mathrm{Re}~\left(\mathrm{InverseFFT}\left(\mathrm{FFTShift}\left(\mathrm{img_fourier_blur}\right)\right)\right):$$\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{Embed}\left(\mathrm{FitIntensity}\left(\mathrm{Create}\left(\mathrm{img_low_pass}\right)\right)\right)$

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.

 >

 >

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).

 > $\mathrm{img}≔\mathrm{Read}\left(\mathrm{cat}\left(\mathrm{kernelopts}\left(\mathrm{datadir}\right),"/images/tree.jpg"\right)\right):\phantom{\rule[-0.0ex]{0.0em}{0.0ex}}\mathrm{Embed}\left(\mathrm{img}\right)$

 > $\mathrm{Embed}\left(\mathrm{ImagePeriodogram}\left(\mathrm{img}\right)\right)$