Accelerating the pace of engineering and science

# Documentation

## Specialized Transforms

### Chirp Z-Transform

The chirp Z-transform (CZT), useful in evaluating the Z-transform along contours other than the unit circle. The chirp Z-transform is also more efficient than the DFT algorithm for the computation of prime-length transforms, and it is useful in computing a subset of the DFT for a sequence. The chirp Z-transform, or CZT, computes the Z-transform along spiral contours in the z-plane for an input sequence. Unlike the DFT, the CZT is not constrained to operate along the unit circle, but can evaluate the Z-transform along contours described by z = AW–ℓ,  = 0, …, M – 1, where A is the complex starting point, W is a complex scalar describing the complex ratio between points on the contour, and M is the length of the transform.

One possible spiral is

```A = 0.8*exp(j*pi/6);
W = 0.995*exp(-j*pi*.05);
M = 91;
z = A*(W.^(-(0:M-1)));
zplane([],z.')```

czt(x,M,W,A) computes the Z-transform of x on these points.

An interesting and useful spiral set is m evenly spaced samples around the unit circle, parameterized by A = 1 and W = exp(-j*pi/M). The Z-transform on this contour is simply the DFT, obtained by

```y = czt(x)
```

czt may be faster than the fft function for computing the DFT of sequences with certain odd lengths, particularly long prime-length sequences.

### Discrete Cosine Transform

The discrete cosine transform (DCT), closely related to the DFT. The DCT's energy compaction properties are useful for applications like signal coding. The toolbox function dct computes the unitary discrete cosine transform, or DCT, for an input vector or matrix. Mathematically, the unitary DCT of an input sequence x is

$y\left(k\right)=w\left(k\right)\sum _{n=1}^{N}x\left(n\right)\mathrm{cos}\frac{\pi \left(2n-1\right)\left(k-1\right)}{N},\text{ }k=1,\dots ,N,$

where

$w\left(k\right)=\left\{\begin{array}{cc}\frac{1}{\sqrt{N}},& k=1,\\ \sqrt{\frac{2}{N},}& 2\le k\le N.\end{array}$

The DCT is closely related to the discrete Fourier transform; the DFT is actually one step in the computation of the DCT for a sequence. The DCT, however, has better energy compaction properties, with just a few of the transform coefficients representing the majority of the energy in the sequence. The energy compaction properties of the DCT make it useful in applications such as data communications.

The function idct computes the inverse DCT for an input sequence, reconstructing a signal from a complete or partial set of DCT coefficients. The inverse discrete cosine transform is

$x\left(n\right)=\sum _{k=1}^{N}w\left(k\right)y\left(k\right)\mathrm{cos}\frac{\pi \left(2n-1\right)\left(k-1\right)}{N},\text{ }n=1,\dots ,N,$

where

$w\left(n\right)=\left\{\begin{array}{cc}\frac{1}{\sqrt{N}},& n=1,\\ \sqrt{\frac{2}{N}},& 2\le n\le N.\end{array}$

Because of the energy compaction mentioned above, it is possible to reconstruct a signal from only a fraction of its DCT coefficients. For example, generate a 25 Hz sinusoidal sequence, sampled at 1000 Hz:

```t = (0:1/999:1);
x = sin(2*pi*25*t);```

Compute the DCT of this sequence and reconstruct the signal using only those components with value greater than 0.1 (64 of the original 1000 DCT coefficients):

```y = dct(x)                   % Compute DCT
y2 = find(abs(y) < 0.9);     % Use 17 coefficients
y(y2) = zeros(size(y2));     % Zero out points < 0.9
z = idct(y);                 % Reconstruct signal w/inverse DCT```

Plot the original and reconstructed sequences:

```subplot(2,1,1)
plot(t,x)
title('Original Signal')
subplot(2,1,2)
plot(t,z)
axis([0 1 -1 1])
title('Reconstructed Signal')```

One measure of the accuracy of the reconstruction is

```norm(x-z)/norm(x)
```

that is, the norm of the difference between the original and reconstructed signals, divided by the norm of the original signal. In this case, the relative error of reconstruction is 0.1443. The reconstructed signal retains approximately 85% of the energy in the original signal.

### Hilbert Transform

The Hilbert transform facilitates the formation of the analytic signal. The analytic signal is useful in the area of communications, particularly in bandpass signal processing. The toolbox function hilbert computes the Hilbert transform for a real input sequence x and returns a complex result of the same length,

```y = hilbert(x)
```

where the real part of y is the original real data and the imaginary part is the actual Hilbert transform. y is sometimes called the analytic signal, in reference to the continuous-time analytic signal. A key property of the discrete-time analytic signal is that its Z-transform is 0 on the lower half of the unit circle. Many applications of the analytic signal are related to this property; for example, the analytic signal is useful in avoiding aliasing effects for bandpass sampling operations. The magnitude of the analytic signal is the complex envelope of the original signal.

The Hilbert transform is related to the actual data by a 90° phase shift; sines become cosines and vice versa. To plot a portion of data (solid line) and its Hilbert transform (dotted line), use

```t = 0:1/1024:1;
x = sin(2*pi*60*t);
y = hilbert(x);

plot(t(1:50),real(y(1:50)))
hold on
plot(t(1:50),imag(y(1:50)))
axis([0 0.05 -1.1 2])
legend('Real Part','Imaginary Part')```

The analytic signal is useful in calculating instantaneous attributes of a time series, the attributes of the series at any point in time. The instantaneous amplitude of the input sequence is the amplitude of the analytic signal. The instantaneous phase angle of the input sequence is the (unwrapped) angle of the analytic signal; the instantaneous frequency is the time rate of change of the instantaneous phase angle. You can calculate the instantaneous frequency using the MATLAB® function, diff.

The Walsh–Hadamard transform is a non-sinusoidal, orthogonal transformation technique that decomposes a signal into a set of basis functions. These basis functions are Walsh functions, which are rectangular or square waves with values of +1 or –1. Walsh–Hadamard transforms are also known as Hadamard (see the hadamard function in the MATLAB software), Walsh, or Walsh-Fourier transforms.

The first eight Walsh functions have these values:

IndexWalsh Function Values
01 1 1 1 1 1 1 1
11 1 1 1 -1 -1 -1 -1
21 1 -1 -1 -1 -1 1 1
31 1 -1 -1 1 1 -1 -1
41 -1 -1 1 1 -1 -1 1
51 -1 -1 1 -1 1 1 -1
61 -1 1 -1 -1 1 -1 1
71 -1 1 -1 1 -1 1 -1

The Walsh–Hadamard transform returns sequency values. Sequency is a more generalized notion of frequency and is defined as one half of the average number of zero-crossings per unit time interval. Each Walsh function has a unique sequency value. You can use the returned sequency values to estimate the signal frequencies in the original signal.

Three different ordering schemes are used to store Walsh functions: sequency, Hadamard, and dyadic. Sequency ordering, which is used in signal processing applications, has the Walsh functions in the order shown in the table above. Hadamard ordering, which is used in controls applications, arranges them as 0, 4, 6, 2, 3, 7, 5, 1. Dyadic or gray code ordering, which is used in mathematics, arranges them as 0, 1, 3, 2, 6, 7, 5, 4.

The Walsh–Hadamard transform is used in a number of applications, such as image processing, speech processing, filtering, and power spectrum analysis. It is very useful for reducing bandwidth storage requirements and spread-spectrum analysis. Like the FFT, the Walsh–Hadamard transform has a fast version, the fast Walsh–Hadamard transform (fwht). Compared to the FFT, the FWHT requires less storage space and is faster to calculate because it uses only real additions and subtractions, while the FFT requires complex values. The FWHT is able to represent signals with sharp discontinuities more accurately using fewer coefficients than the FFT. Both the FWHT and the inverse FWHT (ifwht) are symmetric and thus, use identical calculation processes. The FWHT and IFWHT for a signal x(t) of length N are defined as:

$\begin{array}{l}{y}_{n}=\frac{1}{N}\sum _{i=0}^{N-1}{x}_{i}\mathrm{WAL}\left(n,i\right),\\ {x}_{i}=\sum _{i=0}^{N-1}{y}_{n}\mathrm{WAL}\left(n,i\right),\end{array}$

where i = 0,1, …, N – 1 and WAL(n,i) are Walsh functions. Similar to the Cooley-Tukey algorithm for the FFT, the N elements are decomposed into two sets of N/2 elements, which are then combined using a butterfly structure to form the FWHT. For images, where the input is typically a 2-D signal, the FWHT coefficients are calculated by first evaluating across the rows and then evaluating down the columns.

For the following simple signal, the resulting FWHT shows that x was created using Walsh functions with sequency values of 0, 1, 3, and 6, which are the nonzero indices of the transformed x. The inverse FWHT recreates the original signal.

```x = [4 2 2 0 0 2 -2 0]
y = fwht(x)```
```x =

4     2     2     0     0     2    -2     0

y =

1     1     0     1     0     0     1     0```
`x1 = ifwht(y)`
```x1 =

4     2     2     0     0     2    -2     0```

#### Using Walsh-Hadamard Transform for Spectral Analysis and Compression of ECG Signals

The following example uses an electrocardiogram (ECG) signal to illustrate working with the Walsh-Hadamard transform. ECG signals typically are very large and need to be stored for analysis and retrieval at a future time. Walsh-Hadamard transforms are particularly well-suited to this application because they provide compression and thus require less storage space. They also provide rapid signal reconstruction.

Start with an ECG signal. Replicate it to create a longer signal and insert some additional random noise.

```xe = ecg(512);
xr = repmat(xe,1,8);
x = xr + 0.1.*randn(1,length(xr));
```

Transform the signal using the fast Walsh-Hadamard transform. Plot the original signal and the transformed signal.

```y = fwht(x);

subplot(2,1,1)
plot(x)
xlabel('Sample index')
ylabel('Amplitude')
title('ECG Signal')

subplot(2,1,2)
plot(abs(y))
xlabel('Sequency index')
ylabel('Magnitude')
title('WHT Coefficients')
```

The plot shows that most of the signal energy is in the lower sequency values below approximately 1100. Store only the first 1024 coefficients (out of 4096) and see if the signal can be accurately reconstructed from only these stored coefficients.

```y(1025:length(x)) = 0;
xHat = ifwht(y);

figure
plot(x)
hold on
plot(xHat)
xlabel('Sample index')
ylabel('ECG signal amplitude')
legend('Original','Reconstructed')
```

The reproduced signal is very close to the original but has been compressed to a quarter of the size. Storing more coefficients is a tradeoff between increased resolution and increased noise, while storing fewer coefficients may cause loss of peaks.