X( k + 1) =

kn x(n + 1)W n

1 x(n + 1) = ---N

X(k + 1)Wnkn k=0
In these equations, the series subscripts begin with 1 instead of 0 because of MATLABs vector indexing scheme, and

2 j ------ N

WN = e
Note MATLAB uses a negative j for the fft function. This is an engineering convention; physics and pure mathematics typically use a positive j.
fft, with a single input argument x, computes the DFT of the input vector or matrix. If x is a vector, fft computes the DFT of the vector; if x is a rectangular array, fft computes the DFT of each array column.
For example, create a time vector and signal.
t = (0:1/100:10-1/100); x = sin(2*pi*15*t) + sin(2*pi*40*t); % Time vector % Signal
The DFT of the signal, and the magnitude and phase of the transformed sequence, are then
y = fft(x); % Compute DFT of x m = abs(y); p = unwrap(angle(y)); % Magnitude and phase
To plot the magnitude and phase, type the following commands.
f = (0:length(y)-1)*99/length(y); % Frequency vector plot(f,m); title('Magnitude'); set(gca,'XTick',[60 85]); figure; plot(f,p*180/pi); title('Phase'); set(gca,'XTick',[60 85]);

Magnitude x 10

1.0.5 100
A second argument to fft specifies a number of points n for the transform, representing DFT length.

y = fft(x,n);

In this case, fft pads the input sequence with zeros if it is shorter than n, or truncates the sequence if it is longer than n. If n is not specified, it defaults to the length of the input sequence. Execution time for fft depends on the length, n, of the DFT it performs; see the fft reference page in the MATLAB documentation for details about the algorithm.
The inverse discrete Fourier transform function ifft also accepts an input sequence and, optionally, the number of desired points for the transform. Try the example below; the original sequence x and the reconstructed sequence are identical (within rounding error).
t = (0:1/255:1); x = sin(2*pi*120*t); y = real(ifft(fft(x)));

Overview. 2-2 Filter Requirements and Specification. 2-3 IIR Filter Design. 2-5 Classical IIR Filter Design Using Analog Prototyping. 2-7 Comparison of Classical IIR Filter Types. 2-9 FIR Filter Design. Linear Phase Filters. Windowing Method. Multiband FIR Filter Design with Transition Bands Constrained Least Squares FIR Filter Design. Arbitrary-Response Filter Design. Special Topics in IIR Filter Design Analog Prototype Design. Frequency Transformation. Filter Discretization. Selected Bibliography. 2-17. 2-18. 2-19. 2-23. 2-28. 2-32. 2-38. 2-39. 2-39. 2-42

. 2-46

The Signal Processing Toolbox provides functions that support a range of filter design methodologies. The following sections explain how to apply the filter design tools to Infinite Impulse Response (IIR) and Finite Impulse Response (FIR) filter design problems: Filter Requirements and Specification IIR Filter Design FIR Filter Design Special Topics in IIR Filter Design Selected Bibliography
Filter Requirements and Specification
The goal of filter design is to perform frequency dependent alteration of a data sequence. A possible requirement might be to remove noise above 30 Hz from a data sequence sampled at 100 Hz. A more rigorous specification might call for a specific amount of passband ripple, stopband attenuation, or transition width. A very precise specification could ask to achieve the performance goals with the minimum filter order, or it could call for an arbitrary magnitude shape, or it might require an FIR filter. Filter design methods differ primarily in how performance is specified. For loosely specified requirements, as in the first case above, a Butterworth IIR filter is often sufficient. To design a fifth-order 30 Hz lowpass Butterworth filter and apply it to the data in vector x,
[b,a] = butter(5,30/50); y = filter(b,a,x);
The second input argument to butter specifies the cutoff frequency, normalized to half the sampling frequency (the Nyquist frequency).
Frequency Normalization in the Signal Processing Toolbox All of the filter design functions operate with normalized frequencies, so they do not require the system sampling rate as an extra input argument. This toolbox uses the convention that unit frequency is the Nyquist frequency, defined as half the sampling frequency. The normalized frequency, therefore, is always in the interval 0 f 1. For a system with a 1000 Hz sampling frequency, 300 Hz is 300/500 = 0.6. To convert normalized frequency to angular frequency around the unit circle, multiply by. To convert normalized frequency back to hertz, multiply by half the sample frequency.
More rigorous filter requirements traditionally include passband ripple (Rp, in decibels), stopband attenuation (Rs, in decibels), and transition width (Ws-Wp, in hertz).

Ponesided() d

Spectral Estimation Method Overview
The various methods of spectrum estimation available in the Signal Processing Toolbox can be categorized as follows: Nonparametric methods Parametric methods Subspace methods Nonparametric methods are those in which the estimate of the PSD is made directly from the signal itself. The simplest such method is the periodogram. An improved version of the periodogram is Welchs method [8]. A more modern nonparametric technique is the multitaper method (MTM). Parametric methods are those in which the signal whose PSD we want to estimate is assumed to be output of a linear system driven by white noise. Examples are the Yule-Walker autoregressive (AR) method and the Burg method. These methods estimate the PSD by first estimating the parameters (coefficients) of the linear system that hypothetically generates the signal. They tend to produce better results than classical nonparametric methods when the data length of the available signal is relatively short. Subspace methods, also known as high-resolution methods or super-resolution methods, generate frequency component estimates for a signal based on an eigenanalysis or eigendecomposition of the correlation matrix. Examples are the multiple signal classification (MUSIC) method or the eigenvector (EV) method. These methods are best suited for line spectra that is, spectra of sinusoidal signals and are effective in the detection of sinusoids buried in noise, especially when the signal to noise ratios are low. All three categories of methods are listed in the table below with the corresponding toolbox function names. More information about each function is on the corresponding function reference page. See Parametric Modeling on page 4-12 for details about lpc and other parametric estimation functions.
Functions periodogram pwelch, csd, tfe, cohere pmtm
Periodogram Welch Multitaper)
Power spectral density estimate Averaged periodograms of overlapped, windowed signal sections Spectral estimate from combination of multiple orthogonal windows (or tapers) Autoregressive (AR) spectral estimate of a time-series from its estimated autocorrelation function Autoregressive (AR) spectral estimation of a time-series by minimization of linear prediction errors Autoregressive (AR) spectral estimation of a time-series by minimization of the forward prediction errors Autoregressive (AR) spectral estimation of a time-series by minimization of the forward and backward prediction errors Multiple signal classification Pseudospectrum estimate

1 U = --L

which is independent of the choice of window. The addition of U as a normalization constant ensures that the modified periodogram is asymptotically unbiased.

Welchs Method

An improved estimator of the PSD is the one proposed by Welch [8]. The method consists of dividing the time series data into (possibly overlapping) segments, computing a modified periodogram of each segment, and then averaging the PSD estimates. The result is Welchs PSD estimate. Welchs method is implemented in the Signal Processing Toolbox by the pwelch function. By default, the data is divided into eight segments with 50% overlap between them. A Hamming window is used to compute the modified periodogram of each segment. The averaging of modified periodograms tends to decrease the variance of the estimate relative to a single periodogram estimate of the entire data record. Although overlap between segments tends to introduce redundant information, this effect is diminished by the use of a nonrectangular window, which reduces
the importance or weight given to the end samples of segments (the samples that overlap). However, as mentioned above, the combined use of short data records and nonrectangular windows results in reduced resolution of the estimator. In summary, there is a tradeoff between variance reduction and resolution. One can manipulate the parameters in Welchs method to obtain improved estimates relative to the periodogram, especially when the SNR is low. This is illustrated in the following example. Consider an original signal consisting of 301 samples.
randn('state',1) fs = 1000; % Sampling frequency t = (0:0.3*fs)./fs; % 301 samples A = [2 8]; % Sinusoid amplitudes (row vector) f = [150;140]; % Sinusoid frequencies (column vector) xn = A*sin(2*pi*f*t) + 5*randn(size(t)); periodogram(xn,rectwin(length(xn)),1024,fs);
Periodogram PSD Estimate 10
0 Power Spectral Density (dB/Hz)
We can obtain Welchs spectral estimate for 3 segments with 50% overlap with
pwelch(xn,rectwin(150),75,512,fs);

Welch PSD Estimate 10

5 Power Spectral Density (dB/Hz)
In the periodogram above, noise and the leakage make one of the sinusoids essentially indistinguishable from the artificial peaks. In contrast, although the PSD produced by Welchs method has wider peaks, you can still distinguish the two sinusoids, which stand out from the noise floor. However, if we try to reduce the variance further, the loss of resolution causes one of the sinusoids to be lost altogether.

cohere(xn,yn,256,fs,256,128,'none')

Coherence Function 1

Coherence Function Estimate

0.8 0.6 0.4 0.0

250 Frequency
If the input sequence length nfft, window length window, and the number of overlapping data points in a window numoverlap, are such that cohere operates on only a single record, the function returns all ones. This is because the coherence function for linearly dependent data is one.

Parametric Methods

Parametric methods can yield higher resolutions than nonparametric methods in cases when the signal length is short. These methods use a different approach to spectral estimation; instead of trying to estimate the PSD directly from the data, they model the data as the output of a linear system driven by white noise, and then attempt to estimate the parameters of that linear system. The most commonly used linear system model is the all-pole model, a filter with all of its zeroes at the origin in the z-plane. The output of such a filter for white noise input is an autoregressive (AR) process. For this reason, these methods are sometimes referred to as AR methods of spectral estimation. The AR methods tend to adequately describe spectra of data that is peaky, that is, data whose PSD is large at certain frequencies. The data in many practical applications (such as speech) tends to have peaky spectra so that AR models are often useful. In addition, the AR models lead to a system of linear equations which is relatively simple to solve. The Signal Processing Toolbox offers the following AR methods for spectral estimation: Yule-Walker AR method (autocorrelation method) Burg method
Covariance method Modified covariance method All AR methods yield a PSD estimate given by p 1 P AR ( f ) = --- -----------------------------------------------------------------p 2 fs 2jkf f s 1+ a p ( k )e
The different AR methods estimate the AR parameters ap(k) slightly differently, yielding different PSD estimates. The following table provides a summary of the different AR methods.
Burg Characteristics Covariance Modified Covariance Yule-Walker
Does not apply window to data Minimizes the forward and backward prediction errors in the least squares sense, with the AR coefficients constrained to satisfy the L-D recursion
Does not apply window to data Minimizes the forward prediction error in the least squares sense
Does not apply window to data Minimizes the forward and backward prediction errors in the least squares sense
Applies window to data Minimizes the forward prediction error in the least squares sense (also called Autocorrelation method)

Advantages

High resolution for short data records
Better resolution than Y-W for short data records (more accurate estimates) Able to extract frequencies from data consisting of p or more pure sinusoids

Generalized Cosine Windows
Blackman, Hamming, Hann, and rectangular windows are all special cases of the generalized cosine window. These windows are combinations of sinusoidal sequences with frequencies 0, 2 ( N 1 ) , and 4 ( N 1 ) , where N is the window length. One way to generate them is
ind = (0:n-1)'*2*pi/(n-1); w = A - B*cos(ind) + C*cos(2*ind);
where A, B, and C are constants you define. The concept behind these windows is that by summing the individual terms to form the window, the low frequency peaks in the frequency domain combine in such a way as to decrease sidelobe height. This has the side effect of increasing the mainlobe width. The Hamming and Hann windows are two-term generalized cosine windows, given by A = 0.54, B = 0.46 for Hamming and A = 0.5, B = 0.5 for Hann (C = 0 in both cases). The hamming and hann functions, respectively, compute these windows. Note that the definition of the generalized cosine window shown in the earlier MATLAB code yields zeros at samples 1 and n for A = 0.5 and B = 0.5. The Blackman window is a popular three-term window, given by A = 0.42, B = 0.5, C = 0.08. The blackman function computes this window.

Kaiser Window

The Kaiser window is an approximation to the prolate-spheroidal window, for which the ratio of the mainlobe energy to the sidelobe energy is maximized. For a Kaiser window of a particular length, the parameter controls the sidelobe height. For a given , the sidelobe height is fixed with respect to window length. The statement kaiser(n,beta) computes a length n Kaiser window with parameter beta. Examples of Kaiser windows with length 50 and various values for the beta parameter are
n = 50; w1 = kaiser(n,1); w2 = kaiser(n,4); w3 = kaiser(n,9); [W1,f] = freqz(w1/sum(w1),1,512,2); [W2,f] = freqz(w2/sum(w2),1,512,2); [W3,f] = freqz(w3/sum(w3),1,512,2); plot(f,20*log10(abs([W1 W2 W3]))); grid; legend('beta = 1','beta = 4','beta = 9',3) title('Three Kaiser Window Responses') xlabel('Normalized Frequency (Nyquist = 1)') ylabel('Normalized Magnitude (dB)')
Three Kaiser Window Responses 40 Normalized Magnitude (dB) beta = 1 beta = 4 beta = 9 0.2 0.4 0.6 0.8 Normalized Frequency (Nyquist = 1) 1
As increases, the sidelobe height decreases and the mainlobe width increases. To see how the sidelobe height stays the same for a fixed parameter as the length is varied, try
w1 = kaiser(50,4); w2 = kaiser(20,4); w3 = kaiser(101,4); [W1,f] = freqz(w1/sum(w1),1,512,2); [W2,f] = freqz(w2/sum(w2),1,512,2); [W3,f] = freqz(w3/sum(w3),1,512,2); plot(f,20*log10(abs([W1 W2 W3]))); grid; legend('length = 50','length = 20','length = 101') title('Three Kaiser Window Responses, Beta Fixed') xlabel('Normalized Frequency (Nyquist = 1)') ylabel('Normalized Magnitude (dB)')
Three Kaiser Window Responses, Beta Fixed 20 length = 50 length = 20 length = 101
0 Normalized Magnitude (dB)
0.4 0.6 0.8 Normalized Frequency (Nyquist = 1)

P xy(f) 2 C xy(f) = -----------------------------P xx(f)P yy(f)
x and y must be the same length. nfft specifies the FFT length that cohere uses. This value determines the frequencies at which the coherence is estimated. fs is a scalar that specifies the sampling frequency. window specifies a windowing function and the number of samples cohere uses in its sectioning of the x and y vectors. numoverlap is the number of samples by which the window sections overlap for both x and y. Any arguments that you omit from the end of the parameter list use the default values shown below. Cxy = cohere(x,y) uses the following default values:
nfft = min(256,length(x)) fs = 2 window is a periodic Hann window of length nfft. numoverlap = 0 If x is real, cohere estimates the coherence function at positive frequencies only; in this case, the output Cxy is a column vector of length nfft/2 + 1 for nfft even and (nfft + 1)/2 for n odd. If x or y is complex, cohere estimates the coherence function at both positive and negative frequencies, and Cxy has length nfft.
Cxy = cohere(x,y,nfft) uses the FFT length nfft in estimating the power spectrum for x. Specify nfft as a power of 2 for fastest execution. Cxy = cohere(x,y,nfft,fs) returns a vector f of frequencies at which the function evaluates the coherence. fs is the sampling frequency. f is the same size as Cxy, so plot(f,Cxy) plots the coherence function versus properly scaled frequency. fs has no effect on the output Cxy; it is a frequency scaling multiplier. cohere(x,y,nfft,fs,window) specifies a windowing function and the number of samples per section of the vectors x and y. If you supply a scalar for window, cohere uses a Hann window of that length. The length of the window must be less than or equal to nfft; cohere zero pads the sections if the window length exceeds nfft. cohere(x,y,nfft,fs,window,numoverlap) overlaps the sections of x by numoverlap samples.
Note If you use cohere on two linearly related signals [1] with a single, non-overlapping window, the output for all frequencies is Cxy = 1.
You can use the empty matrix [] to specify the default value for any input argument except x or y. For example,
cohere(x,y,[],[],kaiser(128,5));
uses 256 as the value for nfft and 2 as the value for fs.

[x1,x2] = demod(y,fc,fs,'qam') multiplies y by a cosine and a sine of frequency fc and applies a fifth-order Butterworth lowpass filter using filtfilt. x1 = y.*cos(2*pi*fc*t); x2 = y.*sin(2*pi*fc*t); [b,a] = butter(5,fc*2/fs); x1 = filtfilt(b,a,x1); x2 = filtfilt(b,a,x2);
The default method is 'am'. In all cases except 'ppm' and 'pwm', x is the same size as y. If y is a matrix, demod demodulates its columns.
x = demod(y,fc,fs,'pwm','centered') finds the pulse widths assuming they are centered at the beginning of each period. x is length length(y)*fc/fs.

modulate, vco

dftmtx

7dftmtx

Discrete Fourier transform matrix

A = dftmtx(n)

A discrete Fourier transform matrix is a complex matrix of values around the unit circle, whose matrix product with a vector computes the discrete Fourier transform of the vector.
A = dftmtx(n) returns the n-by-n complex matrix A that, when multiplied into a length n column vector x. y = A*x
computes the discrete Fourier transform of x. The inverse discrete Fourier transform matrix is

Ai = conj(dftmtx(n))/n

In practice, the discrete Fourier transform is computed more efficiently and uses less memory with an FFT algorithm

x = 1:256; y1 = fft(x);

than by using the Fourier transform matrix.
n = length(x); y2 = x*dftmtx(n); norm(y1-y2) ans = 1.8297e-009
dftmtx takes the FFT of the identity matrix to generate the transform matrix. convmtx, fft

7diric

Compute the Dirichlet or periodic sinc function
y = diric(x,n) y = diric(x,n) returns a vector or array y the same size as x. The elements of y are the Dirichlet function of the elements of x. n must be a positive integer.
The Dirichlet function, or periodic sinc function, is
------ ( n 1 ) diric(x, n) = sin ( nx 2 ) --------------------------n sin ( x 2 ) x

x = 0, 2, 4, else

for any nonzero integer n. This function has period 2 for n odd and period 4 for n even. Its peak value is 1, and its minimum value is -1 for n even. The magnitude of this function is (1/n) times the magnitude of the discrete-time Fourier transform of the n-point rectangular window.

fftfilt chooses the key parameters L and nfft in different ways, depending on whether you supply an FFT length n and on the lengths of the filter and signal. If you do not specify a value for n (which determines FFT length), fftfilt chooses these key parameters automatically:
If length(x)is greater than length(b), fftfilt chooses values that minimize the number of blocks times the number of flops per FFT. If length(b) is greater than or equal to length(x), fftfilt uses a single FFT of length
2^nextpow2(length(b) + length(x) - 1)
This essentially computes
y = ifft(fft(B,nfft).*fft(X,nfft))
If you supply a value for n, fftfilt chooses an FFT length, nfft, of 2^nextpow2(n)and a data block length of nfft - length(b) + 1. If n is less than length(b), fftfilt sets n to length(b).

conv, filter, filtfilt

[1] Oppenheim, A.V., and R.W. Schafer, Discrete-Time Signal Processing, Prentice-Hall, 1989.

fftshift

7fftshift
Rearrange the outputs of the FFT functions
y = fftshift(x) y = fftshift(x) rearranges the outputs of fft and fft2 by moving the zero frequency component to the center of the spectral density, which is sometimes a more convenient form.
For vectors, fftshift(x) returns a vector with the left and right halves swapped. For arrays, fftshift(x) swaps quadrants one and three with quadrants two and four. The fftshift function is part of the standard MATLAB language.

For any array X

Y = fft2(x)
has Y(1,1) = sum(sum(X)); the DC component of the signal is in the upper-left corner of the two-dimensional FFT. For

Z = fftshift(Y)

the DC component is near the center of the matrix.

fft, fft2

filter

7filter

Filter data with a recursive (IIR) or nonrecursive (FIR) filter
y = filter(b,a,x) [y,zf] = filter(b,a,x) [.] = filter(b,a,x,zi) [.] = filter(b,a,x,zi,dim) y = filter(b,a,x) filters the data in vector x with the filter described by coefficient vectors a and b to create the filtered data vector y. When x is a matrix, filter operates on the columns of x. When x is an N-dimensional array, filter operates on the first nonsingleton dimension. a(1) cannot be 0, and if a(1) 1, filter normalizes the filter coefficients by a(1). [y,zf] = filter(b,a,x) returns the final values zf of the state vector. When x is a vector, the size of the final condition vector zf is max(length(b),length(a))-1, the number of delays in the filter. When x is a multidimensional array, zf is an s-by-c matrix, where:

[.] = pwelch(x,window,noverlap,.,'range') specifies the range of frequency values. This syntax is useful when x is real. The string 'range' can
be either: 'twosided': Compute the two-sided PSD over the frequency range [0,fs). This is the default for determining the frequency range for complex-valued x. - If you specify fs as the empty vector, [], the frequency range is [0,1). - If you dont specify fs, the frequency range is [0, 2). 'onesided': Compute the one-sided PSD over the frequency ranges specified for real x. This is the default for determining the frequency range for real-valued x. The string 'range' can appear anywhere in the syntax after noverlap.
pwelch(x,.) with no output arguments plots the PSD estimate in dB per
unit frequency in the current figure window.
Estimate the PSD of a signal composed of a sinusoid plus noise, sampled at 1000 Hz. Use 33-sample windows with 32-sample overlap, and the default FFT length, and display the two-sided PSD estimate.
randn('state',0); Fs = 1000; t = 0:1/Fs:.3; x = cos(2*pi*t*200) + randn(size(t)); pwelch(x,33,32,[],Fs,'twosided')
% 200Hz cosine plus noise

Welch PSD Estimate 22

pwelch calculates the power spectral density using Welchs method (see references):
1 The input signal vector x is divided into k overlapping segments according
to window and noverlap (or their default values).
2 The specified (or default) window is applied to each segment of x. 3 An nfft-point FFT is applied to the windowed data. 4 The (modified) periodogram of each windowed segment is computed. 5 The set of modified periodograms is averaged to form the spectrum estimate

S(ej).

6 The resulting spectrum estimate is scaled to compute the power spectral
density as S ( e ) F , where F is: - 2 when you do not supply the sampling frequency - fs when you supply the sampling frequency
The number of segments k that x is divided into is calculated as: Eight if you dont specify window, or if you specify it as the empty vector [] mo k = ------------- if you specify window as a nonempty vector or a scalar lo In this equation, m is the length of the signal vector x, o is the number of overlapping samples (noverlap), and l is the length of each segment (the window length).
pburg, pcov, peig, periodogram, pmcov, pmtm, pmusic, psdplot, pyulear
[1] Hayes, M., Statistical Digital Signal Processing and Modeling, John Wiley & Sons, 1996. [2] Stoica, P., and R.L. Moses, Introduction to Spectral Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1997, pp. 52-54. [3] Welch, P.D, The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms, IEEE Trans. Audio Electroacoustics, Vol. AU-15 (June 1967), pp. 70-73.

g = rc2lar(k) g = rc2lar(k) returns a vector of log area ratio parameters g from a vector of reflection coefficients k. k = [0.3090 0.9801 0.0031 0.0082 -0.0082]; g = rc2lar(k) g = 0.6389 4.6002 0.0062 0.0164 -0.0164

rc2poly

7rc2poly
Convert reflection coefficients to a prediction filter polynomial
a = rc2poly(k) [a,efinal] = rc2poly(k,r0) a = rc2poly(k) converts the reflection coefficients k corresponding to the lattice structure to the prediction filter polynomial a, with a(1) = 1. The output a is row vector of length length(k)+1. [a,efinal] = rc2poly(k,r0) returns the final prediction error efinal based on the zero-lag autocorrelation, r0.
Consider a lattice IIR filter given by reflection coefficients k.
k = [0.3090 0.9800 0.0031 0.0082 -0.0082];
Its equivalent prediction filter representation is given by
a = rc2poly(k) a = 1.0000 0.6148 0.9899 0.0000 0.0032 -0.0082
rc2poly computes output a using Levinsons recursion [1]. The function:
1 Sets the output vector a to the first element of k 2 Loops through the remaining elements of k
For each loop iteration i, a = [a + a(i-1:-1:1)*k(i) k(i)].

3 Implements a = [1 a]

ac2poly, latc2tf, latcfilt, poly2rc, rc2ac, rc2is, rc2lar, tf2latc

7rceps

Real cepstrum and minimum phase reconstruction
y = rceps(x) [y,ym] = rceps(x)
The real cepstrum is the inverse Fourier transform of the real logarithm of the magnitude of the Fourier transform of a sequence.
rceps(x) returns the real cepstrum of the real sequence x. The real cepstrum is a real-valued function. [y,ym] = rceps(x) returns both the real cepstrum y and a minimum phase reconstructed version ym of the input sequence.
rceps is an M-file implementation of algorithm 7.2 in [2], that is, y = real(ifft(log(abs(fft(x)))));
Appropriate windowing in the cepstral domain forms the reconstructed minimum phase signal.
w = [1; 2*ones(n/2-1,1); ones(1 - rem(n,2),1); zeros(n/2-1,1)]; ym = real(ifft(exp(fft(w.*y))));
cceps, fft, hilbert, icceps, unwrap
[1] Oppenheim, A.V., and R.W. Schafer, Digital Signal Processing, Englewood Cliffs, NJ, Prentice-Hall, 1975. [2] Programs for Digital Signal Processing, IEEE Press, New York, 1979.

rectpuls

7rectpuls
Generate a sampled aperiodic rectangle

[1] Rabiner, L.R., and O. Herrmann, The Predictability of Certain Optimum Finite Impulse Response Digital Filters, IEEE Trans. on Circuit Theory, Vol. CT-20, No. 4 (July 1973), pp. 401-408. [2] Rabiner, L.R., and B. Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975, pp. 156-157.

resample

7resample
Change sampling rate by any rational factor
y = resample(x,p,q) y = resample(x,p,q,n) y = resample(x,p,q,n,beta) y = resample(x,p,q,b) [y,b] = resample(x,p,q) y = resample(x,p,q) resamples the sequence in vector x at p/q times the original sampling rate, using a polyphase filter implementation. p and q must be positive integers. The length of y is equal to ceil(length(x)*p/q). If x is a matrix, resample works down the columns of x. resample applies an anti-aliasing (lowpass) FIR filter to x during the resampling process. It designs the filter using firls with a Kaiser window. y = resample(x,p,q,n) uses n terms on either side of the current sample, x(k), to perform the resampling. The length of the FIR filter resample uses is proportional to n; larger values of n provide better accuracy at the expense of more computation time. The default for n is 10. If you let n = 0, resample
performs a nearest-neighbor interpolation
y(k) = x(round((k-1)*q/p)+1)
where y(k) = 0 if the index to x is greater than length(x).
y = resample(x,p,q,n,beta) uses beta as the design parameter for the Kaiser window that resample employs in designing the lowpass filter. The default for beta is 5. y = resample(x,p,q,b) filters x using the vector of filter coefficients b. [y,b] = resample(x,p,q) returns the vector b, which contains the coefficients of the filter applied to x during the resampling process.
Resample a simple linear sequence at 3/2 the original rate.
fs1 = 10; t1 = 0:1/fs1:1; x = t1; y = resample(x,3,2); % % % % Original sampling frequency in Hz Time vector Define a linear sequence Now resample it
t2 = (0:(length(y)-1))*2/(3*fs1); % New time vector plot(t1,x,'*',t2,y,'o',-0.5:0.01:1.5,-0.5:0.01:1.5,':') legend('original','resampled'); xlabel('Time')

1.5 original resampled

0.5 0.5

0.5 Time

Notice that the last few points of the output y are inaccurate. In its filtering process, resample assumes the samples at times before and after the given samples in x are equal to zero. Thus large deviations from zero at the end points of the sequence x can cause inaccuracies in y at its end points. The following two plots illustrate this side effect of resample.
x = [1:10 9:-1:1]; y = resample(x,3,2); subplot(2,1,1); plot(1:19,x,'*',(0:28)*2/3 + 1,y,'o'); title('Edge Effects Not Noticeable'); legend('original','resampled');
x = [10:-1:1 2:10]; y = resample(x,3,2); subplot(2,1,2); plot(1:19,x,'*',(0:28)*2/3 + 1,y,'o') title('Edge Effects Very Noticeable'); legend('original','resampled');

Savitzky-Golay smoothing filters (also called digital smoothing polynomial filters or least squares smoothing filters) are typically used to smooth out a noisy signal whose frequency span (without noise) is large. In this type of application, Savitzky-Golay smoothing filters perform much better than standard averaging FIR filters, which tend to filter out a significant portion of the signals high frequency content along with the noise. Although Savitzky-Golay filters are more effective at preserving the pertinent high frequency components of the signal, they are less successful than standard averaging FIR filters at rejecting noise. Savitzky-Golay filters are optimal in the sense that they minimize the least-squares error in fitting a polynomial to each frame of noisy data.
Use sgolay to smooth a noisy sinusoid and display the result and the first and second derivatives.
N = 4; F = 21; [b,g]=sgolay(N,F); x=5*sin(.4*pi*0:.2:199); y=x+randn(1,996); % Noisy sinusoid for n = (F+1)/2:996-(F+1)/2, % Zero-th order derivative (equivalent to sgolayfilt except % that it doesn't compute transients) z0(n)=g(:,1)'*y(n - (F+1)/2 + 1: n + (F+1)/2 - 1)'; % 1st order derivative z1(n)=g(:,2)'*y(n - (F+1)/2 + 1: n + (F+1)/2 - 1)'; % 2nd order derivative z2(n)=2*g(:,3)'*y(n - (F+1)/2 + 1: n + (F+1)/2 - 1)'; end plot([x(1:length(z0))',y(1:length(z0))',z0']) legend('Noiseless sinusoid','Noisy sinusoid',. 'Smoothed sinusoid') figure plot([diff(x(1:length(z0)+1))',z1']) legend('Noiseless first-order derivative',. 'Smoothed first-order derivative') figure plot([diff(diff(x(1:length(z0)+2)))',z2']) legend('Noiseless second-order derivative',. 'Smoothed second-order derivative')
Note The figures below are zoomed to show more detail.

Zero-th order

Noiseless sinusoid Noisy sinusoid Smoothed sinusoid 6

200 220

First derivative

Second derivative

0.25 0.2 0.15 0.1 0.0.05 0.1 0.15 0.2 0.220 Noiseless secondorder derivative Smoothed secondorder derivative
fir1, firls, filter, sgolayfilt
[1] Orfanidis, S.J., Introduction to Signal Processing, Prentice-Hall, Englewood Cliffs, NJ, 1996.

sgolayfilt

7sgolayfilt

Savitzky-Golay filtering

y = sgolayfilt(x,k,f) y = sgolayfilt(x,k,f,w) y = sgolayfilt(x,k,f,w,dim) y = sgolayfilt(x,k,f) applies a Savitzky-Golay FIR smoothing filter to the data in vector x. If x is a matrix, sgolayfilt operates on each column. The polynomial order k must be less than the frame size, f, which must be odd. If k = f-1, the filter produces no smoothing. y = sgolayfilt(x,k,f,w) specifies a weighting vector w with length f, which contains the real, positive-valued weights to be used during the least-squares minimization. If w is not specified or if it is specified as empty, [], w defaults to an identity matrix. y = sgolayfilt(x,k,f,w,dim) specifies the dimension, dim, along which the filter operates. If dim is not specified, sgolayfilt operates along the first non-singleton dimension; that is, dimension 1 for column vectors and nontrivial matrices, and dimension 2 for row vectors.

References

Algorithm development for the Signal Processing Toolbox has drawn heavily upon the references listed below. All are recommended to the interested reader who needs to know more about signal processing than is covered in this manual.
1 Crochiere, R.E., and L.R. Rabiner. Multi-Rate Signal Processing. Englewood
Cliffs, NJ: Prentice Hall, 1983. Pgs. 88-91.
2 IEEE. Programs for Digital Signal Processing. IEEE Press. New York: John

Wiley & Sons, 1979.

3 Jackson, L.B. Digital Filters and Signal Processing. Third Ed. Boston:
Kluwer Academic Publishers, 1989.
4 Kay, S.M. Modern Spectral Estimation. Englewood Cliffs, NJ: Prentice Hall,
5 Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing.
Englewood Cliffs, NJ: Prentice Hall, 1989.
6 Parks, T.W., and C.S. Burrus. Digital Filter Design. New York: John Wiley

& Sons, 1987.

7 Pratt,W.K. Digital Image Processing. New York: John Wiley & Sons, 1991. 8 Percival, D.B., and A.T. Walden. Spectral Analysis for Physical Applications:
Multitaper and Conventional Univariate Techniques. Cambridge: Cambridge University Press, 1993.
9 Proakis, J.G., and D.G. Manolakis. Digital Signal Processing: Principles,
Algorithms, and Applications. Upper Saddle River, NJ: Prentice Hall, 1996.
1 Rabiner, L.R., and B. Gold. Theory and Application of Digital Signal 0
Processing. Englewood Cliffs, NJ: Prentice Hall, 1975.
1 Welch, P.D. The Use of Fast Fourier Transform for the Estimation of Power 1
Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms. IEEE Trans. Audio Electroacoust. Vol. AU-15 (June 1967). Pgs. 70-73.
Filter Requirements and Specification. 2-2 IIR Filter Design. 2-4 Classical IIR Filter Design Using Analog Prototyping. 2-6 Comparison of Classical IIR Filter Types. 2-8 FIR Filter Design. Linear Phase Filters. Windowing Method. Multiband FIR Filter Design with Transition Bands. Constrained Least Squares FIR Filter Design. Arbitrary-Response Filter Design. Special Topics in IIR Filter Design. Analog Prototype Design. Frequency Transformation. Filter Discretization. 2-16 2-17 2-18 2-22 2-27 2-31 2-37 2-38 2-38 2-41

References. 2-45

Filter Requirements and Specification
The Signal Processing Toolbox provides functions that support a range of filter design methodologies. This chapter explains how to apply the filter design tools to Infinite Impulse Response (IIR) and Finite Impulse Response (FIR) filter design problems. The goal of filter design is to perform frequency dependent alteration of a data sequence. A possible requirement might be to remove noise above 30 Hz from a data sequence sampled at 100 Hz. A more rigorous specification might call for a specific amount of passband ripple, stopband attenuation, or transition width. A very precise specification could ask to achieve the performance goals with the minimum filter order, or it could call for an arbitrary magnitude shape, or it might require an FIR filter. Filter design methods differ primarily in how performance is specified. For loosely specified requirements, as in the first case above, a Butterworth IIR filter is often sufficient. To design a fifth-order 30 Hz lowpass Butterworth filter and apply it to the data in vector x,

[b,a] = butter(5,30/50); y = filter(b,a,x);
The second input argument to butter indicates the cutoff frequency, normalized to half the sampling frequency (the Nyquist frequency).
Frequency Normalization in the Signal Processing Toolbox All of the filter design functions operate with normalized frequencies, so they do not require the system sampling rate as an extra input argument. This toolbox uses the convention that unit frequency is the Nyquist frequency, defined as half the sampling frequency. The normalized frequency, therefore, is always in the interval 0 f 1. For a system with a 1000 Hz sampling frequency, 300 Hz is 300/500 = 0.6. To convert normalized frequency to angular frequency around the unit circle, multiply by. To convert normalized frequency back to Hertz, multiply by half the sample frequency.
More rigorous filter requirements traditionally include passband ripple (Rp, in decibels), stopband attenuation (Rs, in decibels), and transition width (WsWp, in Hertz).

Magnitude

0.4 0.-1 10

10 Frequency(rad/sec)

You can design Butterworth, Chebyshev type I, Chebyshev type II, and elliptic filters that meet this type of performance specification. The toolbox order selection functions estimate the minimum filter order that meets a given set of requirements. To meet specifications with more rigid constraints like linear phase or arbitrary filter shape, use the FIR and direct IIR filter design routines.

IIR Filter Design

The primary advantage of IIR filters over FIR filters is that they typically meet a given set of specifications with a much lower filter order than a corresponding FIR filter. Although IIR filters have nonlinear phase, data processing within MATLAB is commonly performed off-line, that is, the entire data sequence is available prior to filtering. This allows for a noncausal, zero-phase filtering approach (via the filtfilt function), which eliminates the nonlinear phase distortion of an IIR filter. The classical IIR filters, Butterworth, Chebyshev types I and II, elliptic, and Bessel, all approximate the ideal brickwall filter in different ways. This toolbox provides functions to create all these types of classical IIR filters in both the analog and digital domains (except Bessel, for which only the analog case is supported), and in lowpass, highpass, bandpass, and bandstop configurations. For most filter types, you can also find the lowest filter order that fits a given filter specification in terms of passband and stopband attenuation, and transition width(s). The direct filter design function yulewalk finds a filter with magnitude response approximating a desired function. This is one way to create a multiband bandpass filter. You can also use the parametric modeling or system identification functions to design IIR filters. These functions are discussed in the Parametric Modeling section of Chapter 4. The generalized Butterworth design function maxflat is discussed in the section Generalized Butterworth Filter Design on page 2-14. The following table summarizes the various filter methods in the toolbox and lists the functions available to implement these methods.

Kaiser Window

The Kaiser window is an approximation to the prolate-spheroidal window, for which the ratio of the mainlobe energy to the sidelobe energy is maximized. For a Kaiser window of a particular length, the parameter controls the sidelobe height. For a given , the sidelobe height is fixed with respect to window length.
The statement kaiser(n,beta) computes a length n Kaiser window with parameter beta. Examples of Kaiser windows with length 50 and various values for the beta parameter are
n = 50; w1 = kaiser(n,1); w2 = kaiser(n,4); w3 = kaiser(n,9); [W1,f] = freqz(w1/sum(w1),1,512,2); [W2,f] = freqz(w2/sum(w2),1,512,2); [W3,f] = freqz(w3/sum(w3),1,512,2); plot(f,20*log10(abs([W1 W2 W3])))
Three Kaiser Window Responses 0 -20 -40
Normalized Magnitude (dB)
-60 -80 -100 -120 -140 -160 -beta = 1 beta = 2 beta = 3
As increases, the sidelobe height decreases and the mainlobe width increases. To see how the sidelobe height stays the same for a fixed parameter as the length is varied, try
w1 = kaiser(50,4); w2 = kaiser(20,4); w3 = kaiser(101,4); [W1,f] = freqz(w1/sum(w1),1,512,2); [W2,f] = freqz(w2/sum(w2),1,512,2); [W3,f] = freqz(w3/sum(w3),1,512,2); plot(f,20*log10(abs([W1 W2 W3])))
Three Kaiser Window Responses, Beta Fixed 0 length = 50 length = 20 length = 101 -20

-120 0

Kaiser Windows in FIR Design
There are two design formulas that can help you design FIR filters to meet a set of filter specifications using a Kaiser window. To achieve a sidelobe height of dB, the beta parameter is 0.1102 ( 8.7 ), > 50 = 0.5842 ( 21 ) 0.4 + 0.07886 ( 21 ) , 50 21 < 21 0, For a transition width of rad/sec, use the length 8 n = ---------------------- + 1 2.285 Filters designed using these heuristics will meet the specifications approximately, but you should verify this. To design a lowpass filter with cutoff frequency 0.5 rad/sec, transition width 0.2 rad/sec, and 40 dB of attenuation in the stopband, try
[n,wn,beta] = kaiserord([0.4 0.6]*pi,[1 0],[0.01 0.01],2*pi); h = fir1(n,wn,kaiser(n+1,beta),'noscale');
The kaiserord function estimates the filter order, cutoff frequency, and Kaiser window beta parameter needed to meet a given set of frequency domain specifications.
The ripple in the passband is roughly the same as the ripple in the stopband. As you can see from the frequency response, this filter nearly meets the specifications.
[H,f] = freqz(h,1,512,2); plot(f,20*log10(abs(H))), grid
FIR Design using Kaiser Window 20

-20 Magnitude (dB)

-100 0
For details on kaiserord, see the reference description in Chapter 6.

Chebyshev Window

The Chebyshev window minimizes the mainlobe width, given a particular sidelobe height. It is characterized by an equiripple behavior, that is, its sidelobes all have the same height. The chebwin function, with length and sidelobe height parameters, computes a Chebyshev window.

( z z ( 1 ) ) ( z z ( 2 ) )Lz z ( m ) ) ( Z(z) H(z) = ---------- = k -------------------------------------------------------------------------------( z p ( 1 ) ) ( z p ( 2 ) )Lz p ( n ) ) ( P(z) - The Zeros field specifies a variable name or value for the zeros vector z, which contains the locations of m zeros. - The Poles field specifies a variable name or value for the zeros vector p, which contains the locations of n poles. - The Gain field specifies a variable name or value for the gain k. For 2nd Order Sections you specify the filter by its second-order section representation:

H k(z) =

b 0k + b 1k z 1 + b 2k z 2 --------------------------------------------------------1 + a 1k z 1 + a 2k z 2
The SOS Matrix field specifies a variable name or a value for the L-by-6 SOS matrix b 01 b 11 b a 11 a 21 sos = b 02 b 12 b a 12 a 22 M M M b 0L b 1L b 2L M M M 1 a 1L a 2L
whose rows contain the numerator and denominator coefficients bik and aik of the second-order sections of H(z): For every filter, you specify: A variable name or a value for the filters sampling frequency in the Sampling Frequency field
Importing a Spectrum. When you import a spectrum, you specify:
A variable name or a value for the power spectral density (PSD) vector in the PSD field A variable name or a value for the frequency vector in the Freq. Vector field
The PSD values in the PSD vector correspond to the frequencies contained in the Freq. Vector vector; the two vectors must have the same length.
Working with Signals, Filters, and Spectra
When a signal, filter, or spectrum is imported into SPTool or created in SPTool, it is displayed in the appropriate list box, as shown below. Using the Edit menu functions and SPTool buttons, you can edit various properties of the data in SPTool and invoke all of SPTools digital signal processing functions.
Signals list Filters list Spectra list

Signal View button

Spectrum View button

Filter View button

Component Lists in SPTool
Each signal, filter, and spectrum in SPTool is displayed in the appropriate Signals list, Filters list, or Spectra list. Signals are displayed with the signal type [vector] or [array]: - A vector signal ([vector]) has one column of data. - An array signal ([array]) has more than one column of data. Filters are displayed with the filter type [design] or [imported]: - A designed filter ([design]) is a filter that was created using the Filter Designer. This type of filter can also be edited in the Filter Designer. - An imported filter ([imported]) is a filter that was imported from the MATLAB workspace or a file, or one that was edited in the Pole/Zero Editor.
Spectra are displayed with the spectrum type [auto]: An auto-spectrum ([auto]) is a spectrum whose source is a single signal, as opposed to the cross-spectrum of two channels of data. spectrum[auto] is the only spectrum type in SPTool.

The value you enter in the Gain edit box specifies the gain, k. The Coordinates pop-up menu lets you position the poles and zeros using either polar or rectangular measurements. When Polar is selected, you can specify the coordinates of the currently selected pole or zero by entering values in the Mag and Angle edit boxes. The resulting location of the selected pole or zero is z = Me
where M is the magnitude specified by Mag, and is the angle in radians specified by Angle. When Rectangular is selected, you can specify the coordinates of the currently selected pole or zero by entering values in the X and Y edit boxes. The resulting location of the selected pole or zero is z = X + jY The Conjugate pair checkbox controls whether the current selection is a complex conjugate pair of poles or zeros, or a lone pole or zero. When a lone pole or zero is selected on the pole-zero plot, the Conjugate pair checkbox is inactive. If you now enable the Conjugate pair checkbox by clicking it, a new conjugate pole or zero (as appropriate) is created to complete the pair. When a conjugate pair of poles or zeros is selected on the plot (by clicking one of the pair), the Conjugate pair checkbox is active. If you now disable the Conjugate pair checkbox by clicking it, the conjugate pair is broken into two lone poles or zeros at the same positions, and only one of these remains selected. When the Conjugate pair checkbox is selected, the coordinates you specify (in either angular or rectangular form) control the location of the currently selected complex conjugate pair of poles or zeros, with locations z = X jY or z = Me
Designing a New Filter with the Pole/Zero Editor. In general, follow these steps to design a new filter using the Specifications panel parameters:
1 Click-and-drag to select the Pole/Zero Editor
option from the Algorithm pop-up menu.
2 Enter the desired filter gain in the Gain edit box.

Pole Zero

3 Select a pole or zero (or conjugate pair) by clicking

on one of the

Add Poles Tool

Add Zeros Tool

Delete Poles/Zeros Tool

At the bottom of the plot area, there are two additional buttons, Send To Back and Delete All. Use Send To Back when you are unable to select a pole or zero because another pole or zero is in the way. Select the interfering pole or zero on the plot, and press Send To Back to move that pole or zero behind all other symbols on the plot (i.e., last in the stacking order). You can then select the desired pole or zero without the interference of the other symbol.
The Delete All button deletes every pole and zero in the plot area.
NOTE The Delete All operation cannot be undone. Once all the poles and zeros are removed, they cannot be restored.
Designing Finite Impulse Response (FIR) Filters
The Filter Designer provides three options for basic FIR filter design. These options allow you to create FIR filters with standard band configurations (lowpass, highpass, bandpass, or bandstop configurations only). The three options for FIR filter design in the Algorithm pop-up menu are: Equiripple FIR, which accesses the toolbox function remez to create an equiripple FIR filter. Least Squares FIR, which accesses the toolbox function firls to create an FIR filter using the least square design method. Kaiser Window FIR, which accesses the fir1 function to create an FIR filter using a Kaiser window.
Example: FIR Filter Design, Standard Band Configuration
In the following example, use the Kaiser window filter design option:
1 Select Kaiser Window FIR as the filter design from the Algorithm pop-up
2 Select bandpass from the Type pop-up menu as the configuration. 3 Set the filters sampling frequency to 2000 Hz by entering this value in the
Sampling Frequency text box.
4 Press Apply to redraw the response with these settings.
NOTE You must press Apply before you change the following parameters.
5 Check the Minimum Order check box to enable automatic order selection.
6 Set Fp1 to 290 and Fp2 to 525.
These fields respectively define the lower and upper passband edge frequencies in Hertz.
7 Set Fs1 to 200 and Fs2 to 625.
These fields respectively define the lower and upper stopband edge frequencies in Hertz.
8 Set Rp (passband ripple) to 4 and Rs (stopband attenuation) to 30.
Rp and Rs are specified in dB.
9 Press the Apply button.
The Filter Designer calls fir1 to create the filter using a Kaiser window. The Filter Designer updates the magnitude plot to show the new filters magnitude response.

Filter Design Options

When the Minimum Order option is disabled, you can specify parameters that define characteristics unique to certain filter types: For equiripple and least squares filters: the weights for error minimization For Kaiser window filters: the cutoff frequency and parameter of the Kaiser window
Order Selection for FIR Filter Design
As described earlier, the FIR filter design options available through the Filter Designer call the toolbox functions remez, firls, and fir1. In calculating filter order, the Filter Designer uses the same guidelines as the toolbox functions: The Equiripple FIR design option calls the remezord order estimation function to determine a filter order that meets a set of specifications. In some cases, remezord underestimates the filter order n. If the filter does not appear to meet the given specifications using Minimum Order order selection, deselect Minimum Order and manually specify a slightly larger order (n+1 or n+2). The Least Squares FIR design option calls the toolbox function firls. Because the toolbox does not provide an order estimation function for use with firls, you cannot use the Minimum Order option with the Least Squares FIR method. The Kaiser Window FIR design option calls kaiserord, the order estimation function, which sometimes underestimates the filter order n. If the filter does not appear to meet the given specifications using Minimum Order order selection, deselect Minimum Order and manually specify a slightly larger order (n+1 or n+2). All of the FIR filter design options in the Filter Designer require an even filter order for the highpass and bandstop configurations. For more information on order selection with the FIR filter design options, see the reference descriptions of remez, remezord, kaiserord, firls, and fir1 in Chapter 6.
Designing Infinite Impulse Response (IIR) Filters
The Filter Designer lets you design a number of classical IIR filters, including Butterworth, Chebyshev type I, Chebyshev type II, and elliptic filters.
Example: Classical IIR Filter Design
In the following example, design a simple Chebyshev type I filter:
1 Select Chebyshev Type I IIR as the filter design from the Algorithm

pop-up menu.

2 Select highpass from the Type menu as the configuration. 3 Set the filters sample frequency to 2000 Hz by entering this value in the
5 Check the Minimum Order check box. 6 Set Fp (passband edge frequency) to 800 and Fs (stopband edge frequency)

to 700.

Fp and Fs are specified in Hertz.
7 Set Rp (passband ripple) to 2.5 and Rs (stopband attenuation) to 35.
8 Press Apply to draw the magnitude response. 9 Deselect Minimum Order and specify a filter order of 7. Press Apply to

1 Select one or more filters from the Filters list in SPTool. 2 Press View in the Filters panel in SPTool.
The Filter Viewer is activated and the selected filters are loaded into the Filter Viewer and displayed.
Basic Filter Viewer Functions
The Filter Viewer has the following components:
A Plots panel for selecting which subplots display in the main plots window A Rulers panel and line display controls for making signal measurements and comparisons A Frequency Axis panel for specifying x-axis scaling in the main plots window A filter identification panel that displays information about the currently selected filter(s) A main plots (display) area for viewing one or more frequency-domain plots for the selected filter(s) Zoom controls for getting a closer look at filter response characteristics When you first open or activate the Filter Viewer, it displays the default plot configuration for the selected filter(s): View (zoom) controls Line display controls

Filter ID panel

Plots panel, including menus for modifying plot characteristics

Frequency Axis panel

Rulers panel, including controls for measuring filter responses
The filters magnitude and phase plots are displayed. The frequency axis of the plots is set to linear, and the frequency axis range is set to [0,Fs/2].
You can choose to display one or any combination of the six available subplots by using the check boxes in the Plots panel, and you can modify many of the plot display characteristics using the pop-up menus in the Plots panel and the Frequency Axis panel.
File Menu. Use Close from the File menu to close the Filter Viewer. Settings you changed and saved using the Preferences dialog box in SPTool are saved and used the next time you open a Filter Viewer. Window Menu. Use the Window menu to select a currently open MATLAB figure window.
Filter Identification Panel
This panel displays the variable names and the highest sampling frequency of the currently selected filters. To change names or sampling frequencies, use Name or Sampling Frequency from the Edit menu in SPTool.

0.5000

0.7500

1.0000

The bandlimited method uses firls to design an interpolation FIR equivalent to that presented in [1]. The polynomial method uses Lagranges polynomial interpolation formula on equally spaced samples to construct the appropriate filter.

decimate interp resample

Decrease the sampling rate for a sequence (decimation). Increase sampling rate by an integer factor (interpolation). Change sampling rate by any rational factor.
[1] Oetken, Parks, and Schler. New Results in the Design of Digital Interpolators. IEEE Trans. Acoust., Speech, Signal Processing. Vol. ASSP-23 (June 1975). Pgs. 301-309.

6invfreqs

Continuous-time (analog) filter identification from frequency data.
[b,a] [b,a] [b,a] [b,a] [b,a] [b,a] = = = = = = invfreqs(h,w,nb,na) invfreqs(h,w,nb,na,wt) invfreqs(h,w,nb,na,wt,iter) invfreqs(h,w,nb,na,wt,iter,tol) invfreqs(h,w,nb,na,wt,iter,tol,'trace') invfreqs(h,w,'complex',nb,na,.)
invfreqs is the inverse operation of freqs; it finds a continuous-time transfer
function that corresponds to a given complex frequency response. From a laboratory analysis standpoint, invfreqs is useful in converting magnitude and phase data into transfer functions.
[b,a] = invfreqs(h,w,nb,na) returns the real numerator and denominator coefficient vectors b and a of the transfer function
B( s ) b(1)s nb + b(2)s nb 1 + L b(nb + 1) + H(s) = ---------- = --------------------------------------------------------------------------------------------na + a(2)s na 1 + L a(na + 1) A(s) + a(1)s whose complex frequency response is given in vector h at the frequency points specified in vector w. Scalars nb and na specify the desired orders of the numerator and denominator polynomials. Frequency is specified in radians between 0 and , and the length of h must be the same as the length of w. invfreqs uses conj(h) at w to ensure the proper frequency domain symmetry for a real filter.
[b,a] = invfreqs(h,w,nb,na,wt) weights the fit-errors versus frequency. wt is a vector of weighting factors the same length as w. invfreqs(h,w,nb,na,wt,iter) and invfreqs(h,w,nb,na,wt,iter,tol) provide a superior algorithm that guarantees stability of the resulting linear system and searches for the best fit using a numerical, iterative scheme. The iter parameter tells invfreqs to end the iteration when the solution has converged, or after iter iterations, whichever comes first. invfreqs defines convergence as occurring when the

norm of the (modified) gradient vector is less than tol. tol is an optional parameter that defaults to 0.01. To obtain a weight vector of all ones, use
invfreqs(h,w,nb,na,[],iter,tol) invfreqs(h,w,nb,na,wt,iter,tol,'trace') displays a textual progress

report of the iteration.

invfreqs(h,w,'complex',nb,na,.) creates a complex filter. In this case no symmetry is enforced, and the frequency is specified in radians between - and.
When building higher order models using high frequencies, it is important to scale the frequencies, dividing by a factor such as half the highest frequency present in w, so as to obtain well conditioned values of a and b. This corresponds to a rescaling of time. Convert a simple transfer function to frequency response data and then back to the original filter coefficients:
a = [1 4]; b = [3]; [h,w] = freqs(b,a,64); [bb,aa] = invfreqs(h,w,4,5) bb = 1.0000 aa = 1.0000 2.0000 3.0000 2.0000 1.0000 4.0000 2.0000 3.0000 2.0000 3.0000
Notice that bb and aa are equivalent to b and a, respectively. However, aa has poles in the left half-plane and thus the system is unstable. Use invfreqss iterative algorithm to find a stable approximation to the system:
[bbb,aaa] = invfreqs(h,w,4,5,[],30) bbb = 0.6816 aaa = 1.0000 3.4676 7.4060 6.2102 2.5413 0.0001 2.1015 2.6694 0.9113 -0.1218
Suppose you have two vectors, mag and phase, that contain magnitude and phase data gathered in a laboratory, and a third vector w of frequencies. You can convert the data into a continuous-time transfer function using invfreqs:
[b,a] = invfreqs(mag.*exp(j*phase),w,2,3);
By default, invfreqs uses an equation error method to identify the best model from the data. This finds b and a in
by creating a system of linear equations and solving them with MATLABs \ operator. Here A(w(k)) and B(w(k)) are the Fourier transforms of the polynomials a and b, respectively, at the frequency w(k), and n is the number of frequency points (the length of h and w). This algorithm is based on Levi [1]. Several variants have been suggested in the literature, where the weighting function wt gives less attention to high frequencies. The superior (output-error) algorithm uses the damped Gauss-Newton method for iterative search [2], with the output of the first algorithm as the initial estimate. This solves the direct problem of minimizing the weighted sum of the squared error between the actual and the desired frequency response points:

If thresh < 1, there will be no noise eigenvectors. This case is not allowed and gives the following error message:
Noise subspace dimension cannot be zero.
When p < n and thresh 1, p specifies the maximum number of signal eigenvectors. However, the threshold test specified by thresh can also take eigenvectors from the signal subspace and assign them to the noise subspace.
This example analyzes a signal vector xx, assuming that two real signals are present in the signal subspace. In this case, the dimension of the signal subspace is 4 because each real sinusoid is the sum of two complex exponentials:
nn = 0:199; xx = cos(0.257*pi*nn) + sin(0.2*pi*nn) + 0.01*randn(size(nn)); [PP,ff] = pmusic(xx,4);
This example analyzes the same signal vector xx with an eigenvalue cutoff of 10% above the minimum. Setting p = Inf forces the signal/noise subspace
decision to be based on thresh. Use eigenvectors of dimension 7 and a sampling frequency Fs of 8 kHz:
[PP,ff] = pmusic(xx,[Inf,1.1],[],8000,7); % window length = 7
With the third and fourth outputs, by plotting the zeros of the noise-eigenvector polynomials, it is possible to create a Root-MUSIC algorithm, as the following zplane plot illustrates:
[PP,ff,v_noise] = pmusic(xx,4); for kk = 1:size(v_noise,2) rr(:,kk) = roots(v_noise(:,kk)); end zplane(rr)
Assume that RR is a square correlation matrix (for example, 7-by-7):
RR = toeplitz(cos(0.1*pi*[0:6])) + 0.1*eye(7); [PP,ff] = pmusic(RR,4,'corr');
Make an observation matrix xx that is rectangular (100-by-7):
xx = reshape(cos(0.257*pi*(0:699)),7,100) + 0.1*randn(7,100); [PP,ff] = pmusic(xx,4);
Use the same signal, but let pmusic form the 100-by-7 data matrix using its window and overlap inputs. In addition, use a longer FFT:
yy = xx(:); [PP,ff] = pmusic(yy,4,512,[],7,0);
If we set p = 0, all the eigenvectors are assigned to the noise subspace. 'ev' specifies the eigenvector weighting. This turns out to be equivalent to MVDL (Capons MLM):
[PP,ff] = pmusic(RR,0,'ev','corr');
The MUSIC estimate is given by the formula P music(f) = --------------------------------------------------------------- = -----------------------------------------N N H H v k e(f) 2 e H( f ) v k v k e(f) k = p + 1 k = p+1
where N is the dimension of the eigenvectors and vk is the k-th eigenvector of the correlation matrix of the input signal. The integer p is the dimension of the signal subspace, so the eigenvectors vk used in the sum correspond to the smallest eigenvalues and also span the noise subspace. The vector e(f) consists of complex exponentials, so the inner product

H v k e(f)

amounts to a Fourier transform. The second form is preferred for computation because the FFT is computed for each vk and then the squared magnitudes are summed. In the eigenvector method, the summation is weighted by the eigenvalues k of the correlation matrix: 1 P ev(f) = ---------------------------------------------------------N H v k e(f) 2 k k = p + 1

performs a nearest-neighbor interpolation:
y(k) = x(round((k-1)*q/p)+1)
where y(k) = 0 if the index to x is greater than length(x).
y = resample(x,p,q,n,beta) uses beta as the design parameter for the Kaiser window that resample employs in designing the lowpass filter. The default for beta is 5. y = resample(x,p,q,b) filters x with b, a vector of filter coefficients. [y,b] = resample(x,p,q) returns the vector b, which contains the coefficients of the filter applied to x during the resampling process.
Resample a simple linear sequence at 3/2 the original rate:
Fs1 = 10; t1 = 0:1/Fs1:1; x = t1; y = resample(x,3,2); % % % % original sampling frequency in Hz time vector define a linear sequence now resample it
t2 = (0:(length(y)1))*2/(3*Fs1); % new time vector plot(t1,x,'*',t2,y,'o',0.5:0.01:1.5,0.5:0.01:1.5,':') legend('original','resampled'); xlabel('Time')

1.5 original resampled

0.5 0.5

0.5 Time

Notice that the last few points of the output y are inaccurate. In its filtering process, resample assumes the samples at times before and after the given samples in x are equal to zero. Thus large deviations from zero at the end
points of the sequence x can cause inaccuracies in y at its end points. The following two plots illustrate this side effect of resample:
x = [1:10 9:1:1]; y = resample(x,3,2); subplot(2,1,1); plot(1:19,x,'*',(0:28)*2/3 + 1,y,'o'); title('Edge Effects Not Noticeable'); legend('original','resampled'); x = [10:1:1 2:10]; y = resample(x,3,2); subplot(2,1,2); plot(1:19,x,'*',(0:28)*2/3 + 1,y,'o') title('Edge Effects Very Noticeable'); legend('original','resampled');
Edge Effects Not Noticeable Edge Effects Very Noticeable original resampled
5 original resampled 15 20
If p or q are not positive integers, resample gives the appropriate error message:
P must be a positive integer. Q must be a positive integer.
If x is not a vector, resample gives the following error message:
Input X must be a vector.
decimate firls interp interp1 intfilt kaiser spline upfirdn

% Design lowpass filter with cutoff at 1/160th of Fs. hh = fir1(300,2/160); % need a very long lowpass filter yy = upfirdn(xx,hh,160,147);
In this example, the filter design and resampling are separate steps. Note that resample would do both steps as one.
upfirdn uses a polyphase interpolation structure. The number of multiply-add
operations in the polyphase structure is approximately (LhLx-pLx)/q where Lh and Lx are the lengths of h[n] and x[n], respectively. A more accurate flops count is computed in the program, but the actual count is still approximate. For long signals x[n], the formula is quite often exact.
There must be one output argument and at least two input arguments. If either of these conditions are violated, upfirdn gives the appropriate error message:
UPFIRDN needs at least two input arguments. UPFIRDN should have exactly one output argument.
If the arrays are sparse, upfirdn gives the error message
H must be full numeric matrix.
When the input signals are in the columns of a matrix and there are multiple filters also in the columns of a matrix, the number of signals and filters must be the same. If they are not, upfirdn gives the error message
X and H must have the same number of columns, if more than one.
The arguments p and q must be integers. If they are not, upfirdn gives the error message
P and/or Q must be greater than zero
If the arguments p and q are not relatively prime, upfirdn gives the warning message
WARNING (upfirdn) p & q have common factor
conv decimate filter interp intfilt resample
Convolution and polynomial multiplication. Decrease the sampling rate for a sequence (decimation). Filter data with a recursive (IIR) or nonrecursive (FIR) filter. Increase sampling rate by an integer factor (interpolation). Interpolation FIR filter design. Change sampling rate by any rational factor.
[1] Crochiere, R.E., and L.R. Rabiner. Multi-Rate Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 1983. Pgs. 88-91. [2] Crochiere, R.E. A General Program to Perform Sampling Rate Conversion of Data by Rational Ratios. In Programs for Digital Signal Processing. IEEE Press. New York: John Wiley & Sons, 1979. Pgs. 8.2-1 to 8.2-7.
Voltage controlled oscillator.
y = vco(x,Fc,Fs) y = vco(x,[Fmin Fmax],Fs) y = vco(x,Fc,Fs) creates a signal that oscillates at a frequency determined by the real input vector or array x with sampling frequency Fs. Fc is the carrier or reference frequency; when x is 0, y is an Fc Hz cosine with amplitude 1 sampled at Fs Hz. x ranges from -1 to 1, where -1 corresponds to a 0 frequency output, 0 to Fc, and 1 to 2*Fc. y is the same size as x. y = vco(x,[Fmin Fmax],Fs) scales the frequency modulation range so that -1 and 1 values of x yield oscillations of Fmin Hz and Fmax Hz respectively. For best results, Fmin and Fmax should be in the range 0 to Fs/2.