Welch - Maple Help

SignalProcessing

 Welch
 attenuates the effect of noise on a signal and estimates the power spectrum

 Calling Sequence Welch( signal, options )

Parameters

 signal - 1-D rtable or list of data.

Options

 • overlapsize: (optional) Non-negative integer which specifies the target minimum overlap size of the segments. The default is 0.
 • samplerate: (optional) Positive numeric value for the sampling rate. The default is 1.0.
 • segmentsize: (optional) Positive integer for the size of the overlapping segments. The default is the largest power of 2 that is not larger than the size of signal.
 • fftnormalization: (optional) One of none, symmetric, or full, indicates the normalization to be applied when using the Fast Fourier Transform (FFT). The default is symmetric.
 • temperendpoints: (optional) Either true or false, specifies whether the power spectrum is to be tempered at the endpoints. The default is false.
 • window: (optional) Either a list, name, or string, specifies the windowing command to be applied to the overlapping segments. The default is "none" (for no windowing to be applied). If a list is passed, the first element provides the name of the windowing command, and any remaining terms are passed as options to the command.
 • windownormalization (optional) Either true or false, indicates if the windowing function is to be normalized. The default is true.
 • frequencyunit: (optional) Unit which specifies the unit of frequency for the periodogram. The default is Unit(Hz). Either of the forms algebraic or Unit(algebraic) is accepted, and the unit must be convertible to a valid unit of frequency.
 • datapowerscale: (optional) Unit which indicates the scaling, if any, to be applied to the power spectrum Vector. Either of the forms algebraic or Unit(algebraic) is accepted, and the unit must be convertible to a valid unit of power (see below for more details). The default is Unit(1/Hz).
 • plotpowerscale: (optional) Unit which indicates the scaling, if any, to be applied to the power spectrum in the periodogram. Either of the forms algebraic or Unit(algebraic) is accepted, and the unit must be convertible to a valid unit of power (see below for more details). The default is Unit(dB/rad/Hz).
 • periodogramoptions: (optional) Additional plot options to be passed when creating the periodogram. The default is [].
 • output: (optional) The type of output. The supported options are:
 – frequencies: Returns a Vector, of float[8] datatype and length the same as signal, containing the frequencies.
 – power: Returns a Vector of float[8] datatype containing the power spectrum. This is the default.
 – periodogram: Returns a periodogram which displays the power spectrum versus the frequencies.
 – record: Returns a record with the previous options.
 – list of any of the above options: Returns an expression sequence with the corresponding outputs, in the same order.

Description

 • The Welch command takes a 1-D rtable or list signal, and applies Welch's Method to attenuate the effect of noise on the signal, at the expense of frequency resolution, and estimate the power spectrum.
 • The power spectrum is estimated using the following steps:
 1 Divide the signal into segments of equal size with equal or nearly-equal overlaps.
 2 Apply the windowing procedure to each segment.
 3 Take the Discrete Fourier Transform (DFT) of each windowed segment.
 4 Compute the power spectra of each transformed segment.
 5 Compute the estimated overall power spectrum, by averaging the segment power spectra, and scaling the average based on the windowing function and number of segments. Note: The Matrix of all the segment DFTs is known as a Short-Time Fourier Transform (STFT).
 • The values of $a=\mathrm{overlapsize}$, $b=\mathrm{segmentsize}$, and $n=\mathrm{numelems}\left(\mathrm{signal}\right)$ must satisfy $2\le n$, $0\le a$, $2\le b$, $a, and $b\le n$.
 • Since the DFT will be computed for each segment, it is suggested that, for larger signal lengths, segmentsize be a power of 2 and no less than 4, so that the FFT will be utilized.
 • The passed value of overlapsize is used to determine the number $c$ of overlapping segments of size $b=\mathrm{segmentsize}$. The values of $b$ and $c$ determine the smallest possible overlap size $p$ and the excess number $q$ of overlaps of size $p+1$. If we denote $n=\mathrm{numelems}\left(\mathrm{signal}\right)$, then $\mathrm{bc}=n+p\left(c-1-q\right)+\left(p+1\right)q$ for the $1 case. When $c=1$, there are no overlaps and the power spectrum is computed directly from the original signal.
 • The value of window, when not passed as a list, should be the name or string, with or without the Window suffix, that corresponds to the windowing command. For example, to use a Hamming window, you can pass window=Hamming or window="HammingWindow". In both cases, the command SignalProcessing[HammingWindow] will be used internally. Similarly, you can pass window=["Exponential",0.5] or window=[ExponentialWindow,0.5] to use SignalProcessing[ExponentialWindow] with parameter value 0.5.
 • To apply a window to a Vector $V$ of length $n$, the window is first applied to another Vector $W$ of size $n$ and filled with ones, and then $V$ is multiplied element-wise by $W$. When windownormalization=true, $W$ is first normalized with respect to its Root Mean Square (RMS).
 • A power spectrum Vector is determined element-wise from a DFT Vector by squaring the absolute values of each component. When the DFT originated from a real-valued signal, the second half of elements of the power spectrum are discarded (due to symmetry), and the remaining elements of the power spectrum are doubled in value. To scale a power spectrum with the datapowerscale and plotpowerscale options, units which are dimensionally equivalent to the following are accepted:
 • 1: No further scaling is performed.
 • 1/Hz: The power spectrum is divided by $r=\mathrm{samplerate}$.
 • 1/rad/Hz: The power spectrum is divided by $2\mathrm{\pi }r$.
 • dB: Each element $u$ of the power spectrum is replaced with $10{\mathrm{Typesetting}:-\mathrm{_Hold}\left(\left[\mathrm{%log}\right]\right)}_{10}\left(u\right)$.
 • dB/Hz: Each element $u$ of the power spectrum is replaced with $10{\mathrm{Typesetting}:-\mathrm{_Hold}\left(\left[\mathrm{%log}\right]\right)}_{10}\left(\frac{u}{r}\right)$.
 • dB/rad/Hz: Each element $u$ of the power spectrum is replaced with $10{\mathrm{Typesetting}:-\mathrm{_Hold}\left(\left[\mathrm{%log}\right]\right)}_{10}\left(\frac{u}{2\mathrm{\pi }r}\right)$.
 • The frequencies Vector has components defined by ${F}_{i}=\frac{\left(i-1\right)r}{b}$, where $b=\mathrm{segmentsize}$ and $r=\mathrm{samplerate}$. If signal is complex-valued, then both the power spectrum Vector $P$ and $F$ have size $b$. When signal is real-valued, on the other hand, $P$ and $F$ will both have size $m=⌊\frac{b}{2}⌋+1$.
 • The samplerate option can also include a unit of frequency. If a unit is provided, and it differs from frequencyunit, then the sample rate will be converted to use the same unit as frequencyunit.
 • If frequencyunit is different than Unit(Hz), then the labels of the periodogram will use frequencyunit rather than Unit(Hz).
 • When temperendpoints=true, the endpoints of the overall power spectrum, determined by the method described above, are halved.
 • If signal is an rtable of type AudioTools:-Audio, the sample rate is inferred from the attributes. Should samplerate also be passed, it will be overridden.
 • Maple will attempt to coerce the provided signal to a 1-D Vector of either float[8] or complex[8] datatype, and an error will be thrown if this is not possible. For this reason, it is most efficient for the passed input to use this datatype.
 • The input signal cannot have an indexing function, and must use rectangular storage.
 • The Welch command is not thread safe.
 • As the underlying implementation of the SignalProcessing package is a module, it is also possible to use the form SignalProcessing:-Welch to access the command from the package. For more information, see Module Members.

Examples

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

Real-Valued Signal

 • First, define a sampling rate and create Vectors for the times and signal:
 > $n≔8192$
 ${n}{≔}{8192}$ (1)
 > $r≔\mathrm{evalhf}\left(\frac{n-1}{2\mathrm{\pi }}\right)$
 ${r}{≔}{1303.63813886571484}$ (2)
 > $T≔\mathrm{Vector}\left(n,i↦\mathrm{evalhf}\left(\frac{i-1}{r}\right),\mathrm{datatype}=\mathrm{float}\left[8\right]\right)$
 ${{\mathrm{_rtable}}}_{{36893628204687711876}}$ (3)
 > $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)$
 ${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)$ (4)
 > $X≔\mathrm{map}\left[\mathrm{evalhf}\right]\left(g,T\right)$
 ${{\mathrm{_rtable}}}_{{36893628204687700564}}$ (5)
 > $N≔\mathrm{Vector}\left[\mathrm{column}\right]\left(\mathrm{Sample}\left(\mathrm{RandomVariable}\left(\mathrm{Normal}\left(0,2\right)\right),n\right)\right)$
 ${{\mathrm{_rtable}}}_{{36893628204687703444}}$ (6)
 > $Y≔X+N$
 ${{\mathrm{_rtable}}}_{{36893628204687697068}}$ (7)
 • Now, let's compare plots of the signals and power spectra:
 > $\mathrm{plots}:-\mathrm{display}\left(\mathrm{Array}\left(\left[\mathrm{dataplot}\left(T,X,\mathrm{style}=\mathrm{line},\mathrm{view}=\left[0..2\mathrm{\pi },-40..40\right],\mathrm{title}="Original Signal"\right),\mathrm{dataplot}\left(T,Y,\mathrm{style}=\mathrm{line},\mathrm{view}=\left[0..2\mathrm{\pi },-40..40\right],\mathrm{title}="Noisy Signal"\right)\right]\right)\right)$

 > $\mathrm{plots}:-\mathrm{display}\left(\mathrm{Array}\left(\left[\mathrm{Welch}\left(X,\mathrm{samplerate}=r,\mathrm{segmentsize}=n,\mathrm{output}=\mathrm{periodogram},\mathrm{periodogramoptions}=\left[\mathrm{title}="Periodogram of Original Signal"\right]\right),\mathrm{Welch}\left(Y,\mathrm{samplerate}=r,\mathrm{segmentsize}=n,\mathrm{output}=\mathrm{periodogram},\mathrm{periodogramoptions}=\left[\mathrm{title}="Periodogram of Noisy Signal"\right]\right)\right]\right)\right)$

 • Finally, apply Welch's Method with a choice of parameter values, and return the periodogram:
 > $\mathrm{Welch}\left(Y,\mathrm{overlapsize}=512,\mathrm{segmentsize}=1024,\mathrm{window}="Hamming",\mathrm{samplerate}=r,\mathrm{temperendpoints}=\mathrm{true},\mathrm{datapowerscale}=\frac{1}{\mathrm{Hz}},\mathrm{plotpowerscale}=\frac{\mathrm{dB}}{\mathrm{rad}\mathrm{Hz}},\mathrm{output}=\mathrm{periodogram}\right)$

Complex-Valued Signal

 • Consider the following complex-valued signal:
 > $\mathrm{sample_rate}≔5000$
 ${\mathrm{sample_rate}}{≔}{5000}$ (8)
 > $\mathrm{Times}≔\mathrm{Vector}\left(\left[\mathrm{seq}\right]\left(0..1,\frac{1}{\mathrm{sample_rate}}\right),\mathrm{datatype}=\mathrm{float}\left[8\right]\right)$
 ${{\mathrm{_rtable}}}_{{36893628204687664660}}$ (9)
 > $g≔t↦5\cdot \mathrm{cos}\left(1000\cdot \mathrm{\pi }\cdot t\right)-3\cdot I\cdot \mathrm{sin}\left(2000\cdot \mathrm{\pi }\cdot t\right)$
 ${g}{≔}{t}{↦}{5}{\cdot }{\mathrm{cos}}{}\left({1000}{\cdot }{\mathrm{\pi }}{\cdot }{t}\right){-}{3}{\cdot }{I}{\cdot }{\mathrm{sin}}{}\left({2000}{\cdot }{\mathrm{\pi }}{\cdot }{t}\right)$ (10)
 > $\mathrm{Signal}≔\mathrm{Vector}\left(\mathrm{map}\left(g,\mathrm{Times}\right),\mathrm{datatype}=\mathrm{complex}\left[8\right]\right)$
 ${{\mathrm{_rtable}}}_{{36893628204687661044}}$ (11)
 • Here, we will produce a table of possible Welch periodograms using the available scalings:
 > $\mathrm{Scales}≔\left[1,\frac{1}{\mathrm{Hz}},\frac{1}{\mathrm{rad}\mathrm{Hz}},\mathrm{dB},\frac{\mathrm{dB}}{\mathrm{Hz}},\frac{\mathrm{dB}}{\mathrm{rad}\mathrm{Hz}}\right]$
 ${\mathrm{Scales}}{≔}\left[{1}{,}\frac{{1}}{{\mathrm{Hz}}}{,}\frac{{1}}{{\mathrm{rad}}{}{\mathrm{Hz}}}{,}{\mathrm{dB}}{,}\frac{{\mathrm{dB}}}{{\mathrm{Hz}}}{,}\frac{{\mathrm{dB}}}{{\mathrm{rad}}{}{\mathrm{Hz}}}\right]$ (12)
 > $\mathrm{Periodograms}≔\left[\mathrm{seq}\right]\left(\mathrm{Welch}\left(\mathrm{Signal},\mathrm{samplerate}=\mathrm{sample_rate},\mathrm{segmentsize}=64,\mathrm{overlapsize}=32,\mathrm{plotpowerscale}=\mathrm{scale},\mathrm{output}=\mathrm{periodogram}\right),\mathrm{scale}=\mathrm{Scales}\right):$
 > $\mathrm{DocumentTools}:-\mathrm{Tabulate}\left(\left[\mathrm{ListTools}:-\mathrm{LengthSplit}\left(\mathrm{Periodograms},3\right)\right]\right):$

 • The Welch command can also be used to find the power spectrum without attenuation or overlapping segments, by choosing the segment size to be equal to the signal size:
 > $\mathrm{Welch}\left(\mathrm{Signal},\mathrm{samplerate}=\mathrm{sample_rate},\mathrm{segmentsize}=\mathrm{numelems}\left(\mathrm{Signal}\right),\mathrm{output}=\mathrm{power}\right)$
 ${{\mathrm{_rtable}}}_{{36893628204579899260}}$ (13)
 > 

Compatibility

 • The SignalProcessing[Welch] command was introduced in Maple 2021.