Editing Hilbert transform (section)

== Discrete Hilbert transform ==
[[File:Bandpass discrete Hilbert transform filter.tif|thumb|400px|right|''Figure&nbsp;1'': Filter whose frequency response is bandlimited to about 95% of the Nyquist frequency]]
[[File:Highpass discrete Hilbert transform filter.tif|thumb|400px|right|''Figure&nbsp;2'': Hilbert transform filter with a highpass frequency response]]
[[File:DFT approximation to Hilbert filter.png|thumb|400px|right|''Figure&nbsp;3''.]]
[[File:Effect of circular convolution on discrete Hilbert transform.png|thumb|400px|right|''Figure&nbsp;4''. The Hilbert transform of {{math|cos(''ωt'')}} is {{math|sin(''ωt'')}}. This figure shows {{math|sin(ωt)}} and two approximate Hilbert transforms computed by the MATLAB library function, {{mono|hilbert()}}]]
[[File:Discrete Hilbert transforms of a cosine function, using piecewise convolution.svg|thumb|400px|right|''Figure&nbsp;5''. Discrete Hilbert transforms of a cosine function, using piecewise convolution]]
For a discrete function, <math>u[n],</math> with [[discrete-time Fourier transform]] (DTFT), <math>U(\omega)</math>, and '''discrete Hilbert transform''' <math>\widehat u[n],</math> the DTFT of <math>\widehat u[n]</math> in the region {{math|1=−''π'' < ω < ''π''}} is given by''':'''

:<math>\operatorname{DTFT} (\widehat u) = U(\omega)\cdot (-i\cdot \sgn(\omega)).</math>

The inverse DTFT, using the [[Convolution theorem#Functions of a discrete variable (sequences)|convolution theorem]], is''':'''<ref>{{harvnb|Carrick|Jaeger|harris|2011|p=2}}</ref><ref>{{harvnb|Rabiner|Gold|1975|p=71 (Eq 2.195)}}</ref>

:<math>
\begin{align}
\widehat u[n] &= {\scriptstyle \mathrm{DTFT}^{-1}} (U(\omega))\ *\ {\scriptstyle \mathrm{DTFT}^{-1}} (-i\cdot \sgn(\omega))\\
&= u[n]\ *\ \frac{1}{2 \pi}\int_{-\pi}^{\pi} (-i\cdot \sgn(\omega))\cdot e^{i \omega n} \,\mathrm{d}\omega\\
&= u[n]\ *\ \underbrace{\frac{1}{2 \pi}\left[\int_{-\pi}^0 i\cdot e^{i \omega n} \,\mathrm{d}\omega - \int_0^\pi i\cdot e^{i \omega n} \,\mathrm{d}\omega \right]}_{h[n]},
\end{align}
</math>

where

:<math>h[n]\ \triangleq \ 
\begin{cases}
             0, & \text{if }n\text{ even}\\
\frac 2 {\pi n} & \text{if }n\text{ odd}
\end{cases}</math>

which is an infinite impulse response (IIR).

'''Practical considerations'''<ref>{{harvnb|Oppenheim|Schafer|Buck |1999 |p=794-795}}</ref>

'''Method 1:'''  Direct convolution of streaming <math>u[n]</math> data with an FIR approximation of <math>h[n],</math> which we will designate by <math>\tilde h[n].</math>  Examples of truncated <math>h[n]</math> are shown in figures 1 and 2.  '''Fig 1''' has an odd number of anti-symmetric coefficients and is called Type&nbsp;III.<ref>{{harvnb|Isukapalli|p=14}}</ref> This type inherently exhibits responses of zero magnitude at frequencies 0 and Nyquist, resulting in a bandpass filter shape.<ref>{{harvnb|Isukapalli|p=18}}</ref><ref>{{harvnb|Rabiner|Gold|1975|p=172 (Fig 3.74)}}</ref> A Type&nbsp;IV design (even number of anti-symmetric coefficients) is shown in '''Fig 2'''.<ref>{{harvnb|Isukapalli|p=15}}</ref><ref>{{harvnb|Rabiner|Gold|1975|p=173 (Fig 3.75)}}</ref> It has a highpass frequency response.<ref>{{harvnb|Isukapalli|p=18}}</ref> Type III is the usual choice.<ref>{{harvnb|Carrick|Jaeger|harris|2011|p=3}}</ref><ref>{{harvnb|Rabiner|Gold|1975|p=175}}</ref> for these reasons''':'''
* A typical (i.e. properly filtered and sampled) <math>u[n]</math> sequence has no useful components at the Nyquist frequency.
* The Type&nbsp;IV impulse response requires a <math>\tfrac{1}{2}</math> sample shift in the <math>h[n]</math> sequence. That causes the zero-valued coefficients to become non-zero, as seen in ''Figure&nbsp;2''. So a Type&nbsp;III design is potentially twice as efficient as Type&nbsp;IV.
* The group delay of a Type&nbsp;III design is an integer number of samples, which facilitates aligning <math>\widehat u[n]</math> with <math>u[n]</math> to create an [[analytic signal]]. The group delay of Type&nbsp;IV is halfway between two samples.
The abrupt truncation of <math>h[n]</math> creates a rippling (Gibbs effect) of the flat frequency response.  That can be mitigated by use of a window function to taper <math>\tilde h[n]</math> to zero.<ref>{{harvnb|Carrick|Jaeger|harris|2011|p=3}}</ref>

'''Method 2:'''  Piecewise convolution.  It is well known that direct convolution is computationally much more intensive than methods like '''[[Overlap-save method|overlap-save]]''' that give access to the efficiencies of the Fast Fourier transform via the convolution theorem.<ref>{{harvnb|Rabiner|Gold|1975|p=59 (2.163)}}</ref> Specifically, the [[discrete Fourier transform]] (DFT) of a segment of <math>u[n]</math> is multiplied pointwise with a DFT of the <math>\tilde h[n]</math> sequence.  An inverse DFT is done on the product, and the transient artifacts at the leading and trailing edges of the segment are discarded.  Over-lapping input segments prevent gaps in the output stream.  An equivalent time domain description is that segments of length <math>N</math> (an arbitrary parameter) are convolved with the periodic function''':'''                  

:<math>\tilde{h}_N[n]\ \triangleq \sum_{m=-\infty}^\infty \tilde{h}[n - mN].</math>

When the duration of non-zero values of <math>\tilde{h}[n]</math> is <math>M < N,</math> the output sequence includes <math>N-M+1</math> samples of <math>\widehat u.</math>&nbsp; <math>M-1</math> outputs are discarded from each block of <math>N,</math> and the input blocks are overlapped by that amount to prevent gaps.

'''Method 3:'''  Same as method 2, except the DFT of <math>\tilde{h}[n]</math> is replaced by samples of the <math>-i \operatorname{sgn}(\omega)</math> distribution (whose real and imaginary components are all just <math>0</math> or&nbsp;<math>\pm 1.</math>)  That convolves <math>u[n]</math> with a [[periodic summation]]''':'''{{efn-ua
|see {{slink|Convolution_theorem#Periodic_convolution|nopage=y}}, Eq.4b}}

:<math>h_N[n]\ \triangleq \sum_{m=-\infty}^\infty h[n - mN],</math>{{spaces|3}}{{efn-ua
|1=A closed form version of <math>h_N[n]</math> for even values of <math>N</math> is:<ref>{{harvnb|Johansson|p=24}}</ref>
<math display="block">
h_N[n] = \begin{cases}
  \frac{2}{N} \cot(\pi n/N) & \text{for }n\text{ odd},\\
  0 & \text{for }n\text{ even}.
\end{cases}
</math>
}}{{efn-ua
|A closed form version of <math>h_N[n]</math> for odd values of <math>N</math> is''':'''<ref>{{harvnb|Johansson|p=25}}</ref>
<math display="block">h_N[n] = \frac{1}{N} \left(\cot(\pi n/N) - \frac{\cos(\pi n)}{\sin(\pi n/N)}\right),</math>
}}
for some arbitrary parameter, <math>N.</math>  <math>h[n]</math> is not an FIR, so the edge effects extend throughout the entire transform.  Deciding what to delete and the corresponding amount of overlap is an application-dependent design issue.

'''Fig 3''' depicts the difference between methods 2 and 3.  Only half of the antisymmetric impulse response is shown, and only the non-zero coefficients.  The blue graph corresponds to method 2 where <math>h[n]</math> is truncated by a rectangular window function, rather than tapered.  It is generated by a Matlab function, '''hilb(65)'''. Its transient effects are exactly known and readily discarded.  The frequency response, which is determined by the function argument, is the only application-dependent design issue. 

The red graph is <math>h_{512}[n],</math> corresponding to method 3.  It is the inverse DFT of the <math>-i \operatorname{sgn}(\omega)</math> distribution.  Specifically, it is the function that is convolved with a segment of <math>u[n]</math> by the [[MATLAB]] function, '''hilbert(u,512)'''.<ref>{{cite web |author=MathWorks |title=hilbert – Discrete-time analytic signal using Hilbert transform |url=http://www.mathworks.com/help/toolbox/signal/ref/hilbert.html |access-date=2021-05-06 |work=MATLAB Signal Processing Toolbox Documentation}}</ref> The real part of the output sequence is the original input sequence, so that the complex output is an [[Analytic signal|analytic representation]] of <math>u[n].</math> 

When the input is a segment of a pure cosine, the resulting convolution for two different values of <math>N</math> is depicted in '''Fig 4''' (red and blue plots). Edge effects prevent the result from being a pure sine function (green plot). Since <math>h_N[n]</math> is not an FIR sequence, the theoretical extent of the effects is the entire output sequence. But the differences from a sine function diminish with distance from the edges. Parameter <math>N</math> is the output sequence length. If it exceeds the length of the input sequence, the input is modified by appending zero-valued elements. In most cases, that reduces the magnitude of the edge distortions. But their duration is dominated by the inherent rise and fall times of the <math>h[n]</math> impulse response.

'''Fig 5''' is an example of piecewise convolution, using both methods 2 (in blue) and 3 (red dots).  A sine function is created by computing the Discrete Hilbert transform of a cosine function, which was processed in four overlapping segments, and pieced back together. As the FIR result (blue) shows, the distortions apparent in the IIR result (red) are not caused by the difference between <math>h[n]</math> and <math>h_N[n]</math> (green and red in '''Fig 3'''). The fact that <math>h_N[n]</math> is tapered (''windowed'') is actually helpful in this context. The real problem is that it's not windowed enough. Effectively, <math>M=N,</math> whereas the overlap-save method needs <math>M < N.</math>