Digital Signal Processing: Functions and Implementations
Impulse Response
n =- 10:1:30; b = [-1 2 3 6 ...]; A = 1, x = zeros (1,41), x (11) = 1, y1 = filter (B, A, x), stem (n, y1), grid;
Zero-Pole Plot (Zerpol)
function [zeros, poles] = zerpol (B, A); bs = roots (B), as = roots (A); zeros = bs; poles = as; polar (angle (bs), abs (as), 'x') hold off;
Frequency Response (Respfrec)
function [module, phase, frequency] = respfrec (B, A, nfrec) linespace frequency = (0,0.5, nfrec), b = B (length (B): -1:1); a = A (length ( A): -1:1); polyval numerator = (b, exp (-j * 2 * pi * frequency)); polyval denominator = (a, exp (-j * 2 * pi * frequency)); sun = numerator ./ denominator; module = abs (sol) = phase angle (sun); subplor (211), plot (rate widget), grid, xlabel ( 'Frequency') ylabel ( 'Module'), subplot (212) , plot (frequency, unwrap (phase)), grid, xlabel ( 'Frequency') ylabel ( 'Phase (radians)');
FFT Decimation in Time (FFTDT)
function [X, t, f] = fftdt (x), N = length (x) n = 0: N-1, c = clock; flops (0); phasors = exp (-j * 2 * pi * n / N)) if (round (log2 (N)) == log2 (N))
X = mfft (x, phasors)
disp ( 'Time') t = Etim (clock, c), disp ( 'Issue Operations'), f = flop
else
error ( 'No es pot 2');
end;
MFFT Function
function X = (x, phasors), N = length (x) if (N == 2) X ( 1) = x (1) + x (2), X (2) = x (1)-x (2);
else
G = (x (1:2: N), phasors (1:2: N)); mfft H = (x (2:2: N), phasors (1:2: N)), X (1: N / 2) = G + phasors (1: N / 2) .* H; X ((N / 2) +1: N) = G + phasors ((N / 2) +1: N) .* H; end
MDFT Function
function X = (x);
x = x (:), N = length (x) n = [0:1: N-1], for k = 0: N-1
Vk = [exp (-j * 2 * pi * k * n / N)] X (k +1) = Vk * x; FR = [exp (-j * 2 * pi * n / M)] VK1 = FR .* Vk, Vk = VK1; end
Inverse Discrete Fourier Transform (IDFT)
midft sun = function (X) X = conj (X) X = mdft (X) X = conj (X); sun = X / (length (X));
Insert Zeros vs. Interpolation
insertaceros xout = function (xin, L); xout = zeros (1, length (xin, L); xout = zeros (1, length (xin) * L); ind = (1: length (xin)) * L - (L-1); xout (ind) = xin;
The “interp” function, after inserting zeros, inserts a low-pass filter with cutoff frequency fc = 0.5 / L, according to the case.
Remove Samples vs. Decimation
quitamuestras xout = function (xin, M); xout = xin (1: M: end);
Decimate also introduces a filter, which attenuates the pair of tones, before the quitamuestras function. fs >= 2wm (overlap)
Filter vs. Convolution
Filter implements a unidirectional digital filter. y = filter (ba, x) filters the data stored in the vector x with the filter formed by the vectors a and b to create y. If x = impulse, then y = h[n].
Convolution: y = conv (x, h) results in the convolution of vectors x and h. The resulting vector length is length (x) + length (h) - 1. If x and h are vectors of polynomial coefficients, convolving them is equivalent to multiplying the two polynomials.
Convolution is a single FIR filter, while “filter” can be used for both FIR and IIR filters.
Sinusoid Function
function sol = (a, for, fs, fi, t0, tf); 1/fs ts = t = t = t0: ts: tf, n = 1: length (t); sun = A * sin (( 2 * pi * n * fo / fs) + fi);
Delta and Complex Functions
function sol = delta (n); sun = (n == 0);
function sol = complex (fg, fz, A, r), A = 2, R = 0.9, n = 0:20; fg = 90, fz = 45, G = A * exp (j * fg * pi/180) z = r * exp (j * fz * pi/180); sun = G * z ^ n;
Cascade and Parallel Implementations
Cascade 1
x = 0.9 * u[n] a = [1], b1 = [1 / 10 * ones (1,10)], b2 = [1 -1], x1 = filter (b1, a, x) y = filter (b2, a, x1);
Cascade 2
x = 0.9 * u (n), x1 = filter (b2, a, x) y = filter (b1, a, x1),
Parallel
x = 0.9 * u (n), x1 = filter (b1, a, x), x2 = filter (b2, a, x) y = x1 + x2;
Commutative property.
Cascade 2 Equivalent
x = 0.9 * u [n], n =- 10 : 1:100; a = [1], b1 = [1 / 10 * ones (1,10)], h2 = delta (n)-delta (n-1); HC2 = filter (b1, a, h2) b = [0.1 0 0 0 ....- 0.1] y = filter (n, y)