INSTITUTE OF AERONAUTICAL ENGINEERING

DIGITAL SIGNAL AND IMAGE PROCESSING

COURSE CODE
AECD31
BY

K. NAGARAJU
24955A0415

INSTITUTE OF AERONAUTICAL ENGINEERING

DUNDIGAL, HYDERABAD-500 043, TELANGANA,
 DSP & Image Processing — Assignment
Butterworth Filters, IIR Design, Transforms, Wavelets and Image
Restoration

Q1. Butterworth Filter Order for the Given LPF Specifications

Given specifications: passband edge frequency Omega_p = 1 rad/s, stopband edge
frequency Omega_s = 2.33 rad/s, passband ripple A_p = 0.5 dB, stopband
attenuation A_s = 22 dB.

Step 1 — Convert dB specifications to linear power ratios

The Butterworth magnitude-squared response is |H(j*Omega)|^2 = 1 / (1 +
(Omega/Omega_c)^(2N)). The passband condition at Omega = Omega_p requires
the attenuation to be at most A_p dB:

epsilon_p^2 = 10^(A_p/10) - 1 = 10^(0.5/10) - 1

epsilon_p^2 = 10^(0.05) - 1 = 1.1220 - 1 = 0.1220
epsilon_s^2 = 10^(A_s/10) - 1 = 10^(22/10) - 1
epsilon_s^2 = 10^(2.2) - 1 = 158.49 - 1 = 157.49

Step 2 — Compute the filter order N

The order formula for a Butterworth filter derived from the two edge-frequency
conditions is:

N >= log(epsilon_s / epsilon_p) / log(Omega_s / Omega_p)

Numerator: log(sqrt(157.49) / sqrt(0.1220))
= log(12.549 / 0.3493) = log(35.93) = 1.5553
Denominator: log(2.33 / 1) = log(2.33) = 0.3674
N >= 1.5553 / 0.3674 = 4.232

Step 3 — Round up to the nearest integer

Since N must be an integer and must satisfy both specifications, we round up:

N = 5

Verification
With N = 5 and Omega_c selected at the passband edge (Omega_c = Omega_p /
epsilon_p^(1/N) = 1 / (0.1220)^(1/10) = 1.259 rad/s to meet exactly at Omega_p),
the stopband attenuation at Omega_s = 2.33 rad/s is: A_s = 10*log10(1 +
(2.33/1.259)^10) = 10*log10(1 + (1.851)^10) = 10*log10(1 + 370.8) =
10*log10(371.8) = 25.7 dB, which exceeds the required 22 dB. The 5th-order
Butterworth filter satisfies all specifications.

Q2. Conversion of Analog Filter H(s) to Digital IIR Filter Using Impulse
Invariant Method
Given: H(s) = 1 / [(s + 0.2)(s + 0.6)]
Step 1 — Partial Fraction Decomposition
Expand H(s) into partial fractions to identify individual first-order terms (poles):

H(s) = A/(s + 0.2) + B/(s + 0.6)

A = lim_{s -> -0.2} (s+0.2)*H(s) = 1/(-0.2+0.6) = 1/0.4 = 2.5
B = lim_{s -> -0.6} (s+0.6)*H(s) = 1/(-0.6+0.2) = 1/(-0.4) = -2.5
H(s) = 2.5/(s + 0.2) - 2.5/(s + 0.6)

Step 2 — Apply the Impulse Invariant Mapping

In the impulse invariant method, each analog pole at s = -a maps to a digital pole at
z = e^(-aT), where T is the sampling period (assume T = 1 s for normalisation). The
transformation rule for each partial fraction term is:

A/(s + a) --[impulse invariant]--> A / (1 - e^(-aT) * z^(-1))

Step 3 — Compute the Digital Poles

For the first pole: a1 = 0.2, e^(-0.2*1) = e^(-0.2) = 0.8187
For the second pole: a2 = 0.6, e^(-0.6*1) = e^(-0.6) = 0.5488

Step 4 — Write the Digital Transfer Function H(z)

Applying the mapping to each term:

H(z) = 2.5 / (1 - 0.8187*z^(-1)) - 2.5 / (1 - 0.5488*z^(-1))

Step 5 — Combine into a Single Fraction

Let D1 = 1 - 0.8187*z^(-1) and D2 = 1 - 0.5488*z^(-1). Combining:

H(z) = [2.5*D2 - 2.5*D1] / (D1 * D2)

Numerator: 2.5*(1-0.5488*z^(-1)) - 2.5*(1-0.8187*z^(-1))
= 2.5 - 1.372*z^(-1) - 2.5 + 2.047*z^(-1)
= 0.675 * z^(-1)
Denominator: (1 - 0.8187*z^(-1))(1 - 0.5488*z^(-1))
= 1 - 1.3675*z^(-1) + 0.4495*z^(-2)
Final Result: H(z) = 0.675*z^(-1) / (1 - 1.3675*z^(-1) + 0.4495*z^(-2))

Key Notes on Impulse Invariant Method

The impulse invariant method preserves the shape of the impulse response in the
digital domain by ensuring h[n] = T*h_a(nT). It maps the s-plane left-half-plane to the
interior of the unit circle, ensuring stability. The main drawback is aliasing: the digital
frequency response is the sum of periodic replicas of the analog response, causing
overlap (aliasing) unless the analog filter is strictly bandlimited. It is therefore most
suitable for lowpass and bandpass filters but not for highpass or bandstop designs.

Q3. Steps to Design a Chebyshev Type I Filter

The Chebyshev Type I filter is an all-pole IIR filter that achieves a steeper roll-off than
the Butterworth filter of the same order by allowing equiripple behaviour in the
passband while maintaining a monotonically decreasing response in the stopband.

Step 1 — Identify Specifications

Specify passband edge Omega_p, stopband edge Omega_s, maximum passband
ripple A_p (dB), and minimum stopband attenuation A_s (dB). Convert to linear ripple
parameters: epsilon = sqrt(10^(A_p/10) - 1). For example, A_p = 1 dB gives epsilon
= sqrt(10^0.1 - 1) = 0.5088.

Step 2 — Determine the Filter Order N

The order is found from: N >=
cosh^(-1)(sqrt((10^(A_s/10)-1)/(10^(A_p/10)-1))) / cosh^(-1)(Omega_s/Omega_p).
The result is rounded up to the nearest integer to guarantee both specifications are
met simultaneously.

Step 3 — Compute Pole Locations

The N poles of the Chebyshev Type I prototype lowpass filter lie on an ellipse in the s-
plane. Define: phi = (pi/2) + ((2k-1)*pi)/(2N) for k = 1, 2, ..., N; sigma_k = -
sinh(a)*sin(phi_k); Omega_k = cosh(a)*cos(phi_k), where a =
(1/N)*sinh^(-1)(1/epsilon). The complex poles are s_k = sigma_k + j*Omega_k for the
left-half-plane poles only.

Step 4 — Form the Transfer Function

The prototype lowpass transfer function is: H_a(s) = K / product_{k=1}^{floor(N/2)}
(s^2 - 2*sigma_k*s + |s_k|^2) [times (s - s_{(N+1)/2}) for odd N], where K
= Omega_p^N / (2^(N-1)*epsilon) ensures unity DC gain for even N, or adjusts for
odd N. For odd N, there is also a real pole.

Step 5 — Frequency Transformation (if required)

The prototype is a normalised lowpass filter with passband edge at 1 rad/s. To design
a highpass, bandpass, or bandstop Chebyshev filter, apply standard frequency
transformations: Lowpass-to-Highpass: s -> Omega_p/s; Lowpass-to-Bandpass: s ->
(s^2 + Omega_0^2)/(BW*s); Lowpass-to-Bandstop: s -> BW*s/(s^2 + Omega_0^2).

Step 6 — Denormalise (Scale to Actual Frequencies)

Replace s with s/Omega_p to shift the passband edge from 1 rad/s to the actual
Omega_p.

Step 7 — Convert to Digital Domain (if needed)

Apply the bilinear transformation s = (2/T)*(z-1)/(z+1) with frequency pre-warping
(Omega_analog
= (2/T)*tan(omega_digital*T/2)) to obtain the digital IIR Chebyshev filter H(z). The
bilinear transform avoids aliasing entirely, unlike the impulse invariant method.

Properties of the Chebyshev Type I Filter

The Chebyshev polynomial C_N(x) satisfies the recurrence C_N(x) = 2x*C_{N-1}(x) -
C_{N-2}(x) with C_0=1, C_1=x. The magnitude response is |H(j*Omega)|^2 = 1/(1 +
epsilon^2 * C_N^2(Omega/Omega_p)). The response oscillates between 1 and
1/sqrt(1+epsilon^2) in the passband (equiripple), and decreases monotonically in
the stopband. For the same attenuation
requirements, Chebyshev Type I requires a lower order than Butterworth,
resulting in a more efficient design.

Q4. DCT Matrix for N = 4 and Verification of Orthogonality

The Discrete Cosine Transform (DCT-II) matrix C of size N x N has elements defined
as:

C[k,n] = alpha(k) * cos((pi*(2n+1)*k)/(2N)) for k,n = 0,1,...,N-1

where alpha(0) = sqrt(1/N) and alpha(k) = sqrt(2/N) for k >= 1.

Step 1 — Compute the scaling factors for N = 4

alpha(0) = sqrt(1/4) = 0.5
alpha(k) = sqrt(2/4) = sqrt(0.5) = 0.7071 for k=1,2,3

Step 2 — Compute each row of the DCT matrix

Row k=0 (DC component): C[0,n] = 0.5 * cos(0) = 0.5 for all n

C[0,:] = [0.5, 0.5, 0.5, 0.5]

Row k=1:

C[1,n] = 0.7071 * cos(pi*(2n+1)/8) for n=0,1,2,3

C[1,0] = 0.7071*cos(pi/8) = 0.7071*0.9239 = 0.6533
C[1,1] = 0.7071*cos(3pi/8) = 0.7071*0.3827 = 0.2706
C[1,2] = 0.7071*cos(5pi/8) = 0.7071*(-0.3827) = -0.2706
C[1,3] = 0.7071*cos(7pi/8) = 0.7071*(-0.9239) = -0.6533
Row k=2:

C[2,n] = 0.7071 * cos(pi*(2n+1)*2/8) = 0.7071*cos(pi*(2n+1)/4)

C[2,0] = 0.7071*cos(pi/4) = 0.7071*0.7071 = 0.5000
C[2,1] = 0.7071*cos(3pi/4) = 0.7071*(-0.7071) = -0.5000
C[2,2] = 0.7071*cos(5pi/4) = 0.7071*(-0.7071) = -0.5000
C[2,3] = 0.7071*cos(7pi/4) = 0.7071*(0.7071) = 0.5000
Row k=3:

C[3,n] = 0.7071 * cos(pi*(2n+1)*3/8)

C[3,0] = 0.7071*cos(3pi/8) = 0.2706
C[3,1] = 0.7071*cos(9pi/8) = -0.6533
C[3,2] = 0.7071*cos(15pi/8) = 0.6533
C[3,3] = 0.7071*cos(21pi/8) = -0.2706

Step 3 — The Complete 4x4 DCT Matrix

k\n n=0 n=1 n=2 n=3

k=0 0.5000 0.5000 0.5000 0.5000

 k=1 0.6533 0.2706 -0.2706 -0.6533

k=2 0.5000 -0.5000 -0.5000 0.5000

k=3 0.2706 -0.6533 0.6533 -0.2706

Step 4 — Verify Orthogonality: C * C^T = I

For an orthogonal matrix, each row has unit norm and distinct rows are orthogonal.
Check row 0 norm: 0.5^2*4 = 1.0. Check row 1 norm:
0.6533^2+0.2706^2+0.2706^2+0.6533^2 = 0.4268+0.0732+0.0732+0.4268 =
1.0. Check rows 0 and 1 inner product: 0.5*(0.6533+0.2706-0.2706-0.6533) =
0.5*0 = 0. All rows are unit vectors and mutually orthogonal, confirming C*C^T = I.
The DCT matrix is therefore an orthogonal (unitary) matrix, which guarantees energy
preservation (Parseval's theorem) and allows perfect reconstruction via the inverse
DCT = C^T.

Q5. Haar Transform for N = 2 and N = 4 Images

The Haar transform is the simplest wavelet transform. It decomposes a signal or
image into averages and differences at progressively coarser resolutions. It is
lossless (for appropriate scaling) and computationally efficient — requiring only
additions and subtractions.

Haar Basis Matrices

For N = 2, the unnormalised Haar matrix is:

H2 = [ 1 1 ] Normalised: H2 = (1/sqrt(2)) * [ 1 1 ]
[ 1 -1 ] [ 1 -1 ]
For N = 4, the Haar matrix is:

Row Col 0 Col 1 Col 2 Col 3 Description

h0 1 1 1 1 Global average

h1 1 1 -1 -1 Low-freq difference

h2 1 -1 0 0 Left half detail

h3 0 0 1 -1 Right half detail

The normalised N=4 Haar matrix multiplies each row by its appropriate scaling
factor: rows h0 and h1 scale by 1/2, rows h2 and h3 scale by 1/sqrt(2). The scaling
ensures H*H^T = I (orthogonality).

Example: N=2 Haar Transform

Input 1D signal (or image row): f = [1, 3]
F = H2 * f = (1/sqrt(2)) * [1+3, 1-3] = (1/sqrt(2)) * [4, -2]
F = [2*sqrt(2), -sqrt(2)] = [2.828, -1.414]
The first coefficient (2.828) is the scaled average (approximation), and the second (-
1.414) is the scaled difference (detail).
Example: N=4 Haar Transform
Input signal: f = [2, 4, 6, 8]
Step 1 — Apply unnormalised Haar matrix rows then scale:
Row h0: (1/2)*(2+4+6+8) = (1/2)*20 = 10
Row h1: (1/2)*(2+4-6-8) = (1/2)*(-8) = -4
Row h2: (1/sqrt(2))*(2-4+0+0) = (1/sqrt(2))*(-2) = -1.414
Row h3: (1/sqrt(2))*(0+0+6-8) = (1/sqrt(2))*(-2) = -1.414
Haar transform: F = [10, -4, -1.414, -1.414]
The first coefficient (10) captures the global mean intensity. The coefficient -4
captures coarse-scale variation. Coefficients -1.414 and -1.414 capture fine-scale
variations in each half. For a 2D N=4 image, the transform is applied separably: first
along rows, then along columns of the row-transformed matrix, yielding a 4x4 Haar
coefficient matrix with the LL (approximation) subband in the top-left and detail
subbands elsewhere.

Q6. Hadamard Transform — Detailed Explanation and Properties

The Hadamard transform (also called the Walsh-Hadamard transform, WHT) is a
generalisation of the Fourier transform using a matrix whose entries are only +1 and
-1. It requires no multiplications — only additions and subtractions — making it
extremely fast to compute in hardware and software. It is widely used in
communications, signal processing, image processing, and quantum computing.

Definition
The N x N Hadamard matrix H_N (where N = 2^n) is defined recursively by the
Sylvester construction:

H_1 = [1]
H_(2N) = (1/sqrt(2)) * [ H_N H_N ]
[ H_N -H_N ]
For N = 4: H_4 = (1/2) * [[1,1,1,1],[1,-1,1,-1],[1,1,-1,-1],[1,-1,-1,1]]

Forward and Inverse Transform

The forward Hadamard transform of a vector f (length N = 2^n) is: F = H_N * f. Since
H_N is symmetric and orthogonal (H_N * H_N^T = I for the normalised version), the
inverse transform is identical: f = H_N * F. This self-inverse property (H_N^(-1) = H_N
for the unnormalised version with factor 1/N) is a distinctive and computationally
convenient property.

Fast Hadamard Transform (FHT)

Exploiting the recursive structure of the Hadamard matrix, the FHT computes the
transform in O(N log N) operations using only additions and subtractions — no
multiplications. The butterfly structure is analogous to the FFT but simpler. For N = 8,
only 24 additions/subtractions are needed instead of 64 for direct matrix
multiplication.
Properties of the Hadamard Transform
Orthogonality: H_N * H_N^T = N * I (unnormalised) or H_N * H_N^T = I
(normalised). This guarantees perfect reconstruction and energy preservation
(Parseval).

Symmetry: H_N = H_N^T (the matrix is symmetric). Combined with orthogonality,

H_N is both symmetric and orthogonal.

Self-inverse: For the unnormalised transform: (1/N)*H_N*H_N = I, so the inverse uses

the same matrix with a scale factor 1/N.

Real-valued: All transform coefficients are real (even for real inputs), unlike the
DFT which produces complex coefficients.

No multiplications: Only +1 and -1 entries means the transform requires only

additions and subtractions — hardware-friendly.

Energy compaction: Natural images tend to have energy concentrated in low-

sequency coefficients, enabling compression.

Sequency ordering: Unlike frequency ordering in the DFT, Hadamard basis functions
are ordered by sequency (number of sign changes per row), analogous to frequency
content.

Applications
The Hadamard transform is used in spread-spectrum communications (CDMA codes
are Hadamard sequences), image compression (as a fast approximation to the DCT),
error-correcting codes (Reed-Muller codes), and feature extraction. In quantum
computing, the Hadamard gate implements the single-qubit H transform and is
fundamental to creating superposition states.

Q7. Noise Probability Density Functions in Image Processing

Noise in digital images arises from sensor imperfections, quantisation, transmission
errors, and environmental interference. Characterising noise by its probability
density function (PDF) is essential for designing appropriate restoration filters. The
most important noise models in image processing are described below.

1. Gaussian Noise
PDF: p(z) = (1/(sigma*sqrt(2*pi))) * exp(-(z-mu)^2 / (2*sigma^2)). The
Gaussian PDF is bell-shaped, symmetric about the mean mu. Gaussian noise arises
from electronic thermal noise, amplifier noise, and readout noise in sensors. It is the
most common noise model because, by the Central Limit Theorem, the sum of many
independent noise sources is approximately Gaussian. It affects all pixel intensities
uniformly. The parameter sigma controls the spread (standard deviation); larger
sigma means more severe noise.

2. Rayleigh Noise
PDF: p(z) = (2/b)*(z-a)*exp(-(z-a)^2/b) for z>=a, 0 otherwise. The Rayleigh
distribution is skewed right (positive skew). It commonly occurs in range imaging
(radar, sonar) and phase images. Mean = a + sqrt(pi*b/4); Variance = b*(4-pi)/4.
3. Erlang (Gamma) Noise
PDF: p(z) = (a^b * z^(b-1) * exp(-az)) / (b-1)! for z>=0. The Erlang distribution
generalises the exponential distribution. It is used to model noise in laser imaging
and certain quantum detectors. Mean = b/a; Variance = b/a^2. For b=1 it reduces to
the exponential distribution.

4. Exponential Noise
PDF: p(z) = a*exp(-az) for z>=0, 0 otherwise. A special case of the Erlang/Gamma
distribution (b=1). Mean = 1/a; Variance = 1/a^2. It models certain types of laser
speckle noise.

5. Uniform Noise
PDF: p(z) = 1/(b-a) for a<=z<=b, 0 otherwise. Uniform noise is flat within [a, b], with
all intensity values equally probable. Mean = (a+b)/2; Variance = (b-a)^2/12. It
arises from quantisation errors and is the best-case noise model in the sense of
maximum entropy for bounded distributions.

6. Impulse (Salt and Pepper) Noise

PDF: p(z) = Pa for z=a (pepper/dark), Pb for z=b (salt/bright), and 1-Pa-Pb for
z=other. Also called bipolar impulse noise or shot noise. It manifests as randomly
occurring white and black pixels. Caused by bit errors in transmission, faulty memory
cells, or ADC errors. Unlike other noise types, its PDF is not a smooth function but
consists of impulses at the extreme intensity values.

PDF Shape Summary Table

Noise Type PDF Shape Key Parameter Typical Source

Gaussian Bell curve, Mean mu, Std Thermal / electronic noise

symmetric sigma
Rayleigh Right-skewed a (offset), b Radar / range imaging
unimodal (spread)
Erlang Right-skewed, multi- l a, b (shape) Laser imaging
moda
Exponential Monotone Rate a Speckle noise
decreasing
Uniform Flat rectangle Bounds a, b Quantisation error

Impulse Two spikes at Probabilities Pa, Transmission errors

extremes Pb

Q8. 2D Continuous Wavelet Transform (CWT) and Comparison with DWT

The Continuous Wavelet Transform (CWT) decomposes a signal or image into
components at all possible scales and translations simultaneously, providing a time-
frequency (or space-frequency) representation with both spatial and frequency
localisation. Unlike the Fourier transform, which uses infinite sinusoids, the CWT uses
localised, oscillatory functions called wavelets.

1D CWT — Foundation
The 1D CWT of a signal f(t) with respect to a mother wavelet psi(t) is: CWT_f(a, b) =
(1/sqrt(|a|)) * integral f(t) * psi*((t-b)/a) dt, where a is the scale parameter (a > 0,
controls frequency: small a = high frequency), and b is the translation parameter
(controls spatial position). The factor 1/sqrt(|a|) preserves energy across scales.
2D Continuous Wavelet Transform
The 2D CWT extends the 1D concept to images f(x, y). The 2D mother wavelet psi(x,
y) must be square-integrable. The 2D CWT is:

CWT_f(a, b_x, b_y, theta) = (1/a) * double_integral f(x,y) *

psi*((r_theta(x-b_x, y-b_y))/a) dx dy

where a > 0 is scale, (b_x, b_y) is the 2D translation vector, and theta is the
orientation angle (rotation). The rotation operator r_theta rotates the wavelet so that
it can detect edges and features at any orientation in the image. This gives the 2D
CWT four degrees of freedom: scale a, two translations b_x, b_y, and rotation theta.

The 2D CWT with isotropic wavelets (such as the Mexican Hat or Morlet wavelet)
detects blob-like features at various scales. Directional wavelets (elongated, oriented
wavelets) detect oriented edges, ridges, and textures. Common 2D wavelets include
the 2D Morlet wavelet (modulated Gaussian), 2D Mexican Hat (Laplacian of
Gaussian), and curvelets/shearlets for anisotropic analysis.

Comparison: CWT vs DWT

Property 2D CWT 2D DWT

Scale sampling Continuous: all a > 0 Dyadic: a = 2^j, j integer

Translation Continuous: all (b_x, b_y) Discrete lattice points

Orientation Continuous: all theta Fixed: H, V, Diagonal (3

subbands)
Redundancy Highly redundant Critically sampled (non-
redundant)
Computation Very expensive O(N) fast filter bank
algorithm
Reconstruction Via admissibility condition act with perfect reconstruction
Exfilters
Edge detection Excellent at all angles Limited to 3 orientations

Shift invariance Yes (due to redundancy) No (shift-variant)

Memory Very large Compact

Applications Feature detection, texture Compression (JPEG2000),

analysis denoising
Coefficients Complex (analytic Real or complex
wavelet)

In practice, the DWT is overwhelmingly used in image compression and denoising

due to its computational efficiency, while the CWT is used in analysis tasks where
complete representation across scales and orientations is more important than
efficiency.

Q9. Power Density Functions of Uniform, Salt and Pepper, and Gaussian
Noise
The power spectral density (PSD) and probability density function (PDF) of noise
characterise its statistical and spectral behaviour. Here we describe the PDF forms
and their characteristic shapes for the three most important noise types in image
restoration.
1. Uniform Noise PDF
p(z) = 1/(b-a) for a <= z <= b; 0 otherwise.
The PDF is a perfect rectangle of height 1/(b-a) between intensities a and b. Mean =
(a+b)/2. Variance = (b-a)^2/12. The power spectrum of uniform noise is flat (white
noise) — all frequency components have equal power. This means the noise is
spectrally unstructured and affects all spatial frequencies equally in an image. On a
plot, the PDF appears as a flat-topped rectangle, with the height inversely
proportional to the range (b-a). Wider range means lower but wider PDF.

2. Salt and Pepper (Impulse) Noise PDF

p(z) = Pa if z = a (pepper, dark), Pb if z = b (salt, bright), 1-Pa-Pb
for intermediate z.

The PDF consists of two Dirac delta impulses (mathematically) at the extreme
intensity values — one at near-zero (pepper: black pixels) with probability Pa, and
one at near-maximum intensity (salt: white pixels) with probability Pb. All other
intensities have much smaller probability (or zero in pure impulse noise). On a plot,
this appears as two tall vertical spikes at the extreme ends of the intensity axis. If Pa
= Pb, the noise is symmetric. Salt and pepper noise has a flat power spectrum as
well, since the corrupted pixels occur at random spatial locations.

3. Gaussian Noise PDF

p(z) = (1/(sigma*sqrt(2*pi))) * exp(-(z-mu)^2/(2*sigma^2))
The Gaussian PDF is a smooth bell curve centred at the mean mu with spread
controlled by the standard deviation sigma. About 68% of noise values fall within one
sigma of the mean, 95% within two sigma, and 99.7% within three sigma. The PDF is
symmetric about mu. On a plot, it appears as the classic Gaussian bell shape —
tallest at z = mu, falling off rapidly for |z - mu| >> sigma. Narrower curves indicate
lower noise variance; broader curves indicate higher noise. Gaussian noise is the
most theoretically tractable and commonly assumed noise model.

PDF Shape Sketch Summary

Uniform: Rectangle between [a, b], flat top, zero elsewhere. Gaussian: Bell curve,
peak at z=mu, symmetric, extends to +/-infinity. Salt and Pepper: Two sharp vertical
spikes at z=a (pepper) and z=b (salt), near-zero elsewhere.

Noise PDF Shape Description Mean Variance Spectrum

Uniform Flat rectangle [a,b] (a+b)/2 (b-a)^2/12 White (flat)

Gaussian Bell curve at mu mu sigma^2 White (flat)

Salt & Pepper Two impulses at a, b a+Pb*b)/ bH) igh (bimodal) White (flat)
(Pa* (Pa+P

Q10. Mean Filter and Geometric Mean Filter for Image Restoration
Spatial filtering is one of the most fundamental approaches to image restoration.
Order-statistic and averaging filters work in a local neighbourhood (kernel) around
each pixel to estimate the clean pixel value, effectively suppressing noise by
exploiting the spatial correlation of natural
images.

A. Arithmetic Mean (Average) Filter

The arithmetic mean filter replaces each pixel f_hat(x,y) with the arithmetic
average of all pixel values in the rectangular kernel S_xy of size m x n:

f_hat(x,y) = (1/(m*n)) * sum_{(s,t) in S_xy} g(s,t)

where g(s,t) is the noisy observed image. This filter is the simplest and most intuitive
noise reduction method. It is optimal (in the mean-squared error sense) for removing
additive Gaussian white noise — the averaging operation exploits the zero-mean
property of Gaussian noise, so over many samples the noise contributions cancel
while the signal average is preserved.

Properties of the Arithmetic Mean Filter

Noise reduction: Reduces zero-mean additive noise — the standard deviation of the
filtered noise is sigma/sqrt(m*n), so a 3x3 kernel reduces noise standard deviation by
a factor of 3. Blurring: Larger kernels reduce more noise but cause more blurring
(loss of fine detail, edge sharpening, and texture). Linear operation: The arithmetic
mean filter is a linear shift-invariant (LSI) filter, equivalent to convolution with an m x
n box kernel — computationally very efficient. Effect on impulse noise: The arithmetic
mean filter is poor for impulse (salt and pepper) noise because a single extreme-
valued pixel in the neighbourhood strongly pulls the average away from the correct
value. Boundary handling: At image borders, zero-padding, replicate-padding, or
cropping are used to handle incomplete neighbourhoods.

B. Geometric Mean Filter

The geometric mean filter computes the geometric mean (product raised to the
power 1/(m*n)) of pixel values in the neighbourhood S_xy:

f_hat(x,y) = [ product_{(s,t) in S_xy} g(s,t) ]^(1/(m*n))

In logarithmic form: f_hat(x,y) = exp( (1/(m*n)) * sum_{(s,t)}
ln(g(s,t)) )

Properties of the Geometric Mean Filter

Noise reduction: The geometric mean filter also smooths noise but tends to lose less
image detail compared to the arithmetic mean filter of the same neighbourhood size.
This is because the geometric mean is always less than or equal to the arithmetic
mean (AM-GM inequality), with equality only when all values in the neighbourhood
are identical. For natural images with gradual spatial variation, this results in less
blurring while still averaging out noise.

Nonlinearity: Unlike the arithmetic mean filter, the geometric mean filter is a non-
linear operation (due to the product and power). This non-linearity gives it better
performance for certain noise distributions, particularly multiplicative noise (speckle
noise), where noise is modelled as f * eta (product of signal and noise). Taking the
logarithm converts multiplicative noise to additive noise, and the geometric mean
filter can then suppress it while the arithmetic mean cannot.

Handling zero values: The geometric mean filter fails when any pixel in the
neighbourhood is zero, since the product becomes zero regardless of other values.
A common remedy is to add a
small constant epsilon before filtering: replace g(s,t) with g(s,t) + epsilon.

Comparison summary: The arithmetic mean filter is optimal for Gaussian (additive)
noise and is simpler to implement; it tends to blur edges more. The geometric mean
filter is better for multiplicative noise and log-normal distributions, preserves edges
better at the cost of slightly higher computational complexity (requiring logarithms
and exponentials).

Property Arithmetic Mean Geometric Mean

Formula Average of g(s,t) (Product g(s,t))^(1/mn)

Noise model Additive Gaussian Multiplicative / Log-normal

Linearity Linear (LSI) Nonlinear

Blurring More blurring Less blurring

Zero handling No problem Fails if any pixel = 0

Computation Simple (additions) Log, multiply, exp

Impulse noise Poor Poor (but slightly better)

Implementation Box filter (fast) Logarithm domain

— End of Assignment —