Fourier Analysis and Filtering
|
This content was generated with the assistance of AI and should be verified against the official documentation before being relied on in production. The reference material consulted while preparing these pages is listed in the section bibliography. |
Written against MATLAB R2026b. The snippets in the "Octave 11.3.0 output" blocks were run in GNU Octave 11.3.0 using only code that Octave also accepts (fft, ifft, fft2, fftshift, ifftshift, conv, deconv, filter); the MATLAB layout of the same numbers may differ in spacing. Everything else is labelled "not run here" and was checked line by line against the R2026b documentation. The signals are built from sines, not random numbers, so the results are deterministic in both products.
The Fourier transform decomposes data, seen as a function of time or space, into frequency components. In MATLAB the discrete version is computed by fft, a fast algorithm whose cost is much lower than the direct sum. The same documentation family covers convolution (conv) and digital filtering (filter), which modify the amplitude or phase of data through a transfer function. Base MATLAB has all of these; spectral estimation, filter design and zero-phase filtering belong to the Signal Processing Toolbox (see The boundary with Signal Processing Toolbox). Smoothing without any toolbox is also covered by movmean and smoothdata, described with the data-preprocessing functions on Descriptive Statistics and Preprocessing.
The heatflow signal used on this page
A thermocouple on the plate samples the temperature at 10 Hz (fs = 10) for 10 s, so N = 100 samples. The synthetic signal has a 20 degree mean, a 1 Hz heater cycle of amplitude 2 and a 3 Hz vibration of amplitude 0.5. The scenario is described under The heatflow scenario; the time column of sensors.csv is a datetime, so a real log is first turned into seconds with seconds(t - t(1)) and must be evenly spaced (otherwise use retime on a timetable, see Timetables).
fft, ifft and their conventions
Y = fft(X)
Y = fft(X,n)
Y = fft(X,n,dim)fft(X) returns the discrete Fourier transform, the same size as X. For a vector it transforms the vector; for a matrix it transforms each column; for an array it works along the first dimension whose size is not 1. fft(X,n) pads with trailing zeros, or truncates, to n points. fft(X,n,dim) works along dimension dim, so fft(X,n,2) transforms each row. The transform is the sum
with no scaling factor, so a constant signal of value c and length N gives Y(1) = N*c. ifft computes the inverse, including the division by N, so ifft(fft(x)) returns x up to rounding error. The result of fft is complex even for real input. Real input produces a conjugate-symmetric spectrum: the second half repeats the first half as complex conjugates, which is why the single-sided spectrum below uses only the first half. When ifft receives a spectrum that is symmetric up to rounding, it can return a tiny imaginary part; the 'symmetric' option of ifft removes it (see the reference page, not run here).
Source: MATLAB documentation: fft
Frequency vector and the single-sided amplitude spectrum
Element k (counting from 1) of fft(x) belongs to the frequency (k-1)*fs/N. The first element is the zero frequency (the mean times N), the first half holds the positive frequencies up to the Nyquist frequency fs/2, and the second half the negative ones. The resolution is fs/N. The documentation’s recipe to read amplitudes divides by N, takes the absolute value, keeps the first N/2+1 points and doubles all but the first and last, which have no conjugate partner:
fs = 10; N = 100;
t = (0:N-1)'/fs;
x = 20 + 2*sin(2*pi*1*t) + 0.5*sin(2*pi*3*t);
Y = fft(x);
P2 = abs(Y/N);
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:N/2)/N;
fprintf('%5.1f Hz %6.3f\n', [f(P1>1e-6); P1(P1>1e-6)']);
disp(size(Y)); disp(isreal(Y)); disp(max(abs(x - ifft(Y))) < 1e-12)Octave 11.3.0 output:
0.0 Hz 20.000
1.0 Hz 2.000
3.0 Hz 0.500
100 1
0
1The three components are recovered with their exact amplitudes, Y is a 100-by-1 complex column, and the round trip is exact to rounding. This assumes an even N; for an odd length the Nyquist element does not exist and the slice is P2(1:(N+1)/2).
fftshift and the negative frequencies
Y = fftshift(X) swaps the left and right halves of a vector (the first and third quadrants of a matrix), moving the zero frequency to the centre; fftshift(X,dim) works along one dimension. With an odd length the middle element counts as part of the left half. ifftshift undoes it. For an even length N the shifted frequency axis is fs*(-N/2:N/2-1)/N.
fs = 10; N = 8;
f = fs*(-N/2:N/2-1)/N; disp(f)
disp(fftshift(0:7)); disp(ifftshift(fftshift(0:7)))Octave 11.3.0 output:
-5.0000 -3.7500 -2.5000 -1.2500 0 1.2500 2.5000 3.7500
4 5 6 7 0 1 2 3
0 1 2 3 4 5 6 7Source: MATLAB documentation: fftshift
Spectral leakage and zero padding
The transform assumes that the N samples are one period of a periodic signal. If the record holds a whole number of periods of every component, each component lands on one bin. Otherwise its energy spreads over many bins ("leakage"). With the sampling above, a 1 Hz sine fits ten whole periods in 10 s, a 1.25 Hz sine does not:
N = 100; fs = 10; t = (0:N-1)'/fs;
x = sin(2*pi*1.25*t);
P = abs(fft(x)/N); disp(sum(P(1:N/2+1) > 0.01))
x = sin(2*pi*1*t);
P = abs(fft(x)/N); disp(sum(P(1:N/2+1) > 0.01))Octave 11.3.0 output:
25
1The off-bin sine puts more than 0.01 of amplitude in 25 bins, the on-bin one in a single bin. Windowing reduces leakage (hann, hamming and friends are Signal Processing Toolbox functions). Zero padding with fft(x,n) does not add information or improve the resolution fs/N of the original record; it interpolates the spectrum on a finer grid and lets you choose a fast length. The documentation’s Gaussian pulse example pads a record of 44,101 points (a large prime) to 2^nextpow2(L), because the execution time of fft depends on the length: lengths with only small prime factors are much faster than lengths with large prime factors.
Two dimensions: fft2
Y = fft2(X) is the 2-D transform of a matrix, equivalent to fft(fft(X).').'. Written with the transpose operator that way, note the .' is the non-conjugating transpose. fft2(X,m,n) pads or truncates to m-by-n. For an array with more than two dimensions it transforms each 2-D page. fftshift centres the zero frequency in the 2-D result, which is how the documentation displays images and repeated blocks with imagesc(abs(fftshift(Y))). An impulse in the plate temperature field has a flat spectrum, a quick check of the transform:
A = zeros(4,4); A(2,2) = 1; B = fft2(A);
disp(abs(B))
disp(max(max(abs(ifft2(B) - A))) < 1e-12)Octave 11.3.0 output:
1 1 1 1
1 1 1 1
1 1 1 1
1 1 1 1
1For the plate’s temperature field T (see Matrices and Arrays), abs(fftshift(fft2(T - mean(T,"all")))) shows which spatial wavelengths dominate (not run here). fftn and ifftn extend this to N dimensions (reference pages not opened for this page).
Source: MATLAB documentation: fft2
Convolution with conv
w = conv(u,v)
w = conv(u,v,shape)conv returns the convolution of two vectors, of length length(u)+length(v)-1. The documentation points out that, algebraically, convolution is polynomial multiplication: the convolution of the coefficients of and is [2 7 2 7], the coefficients of . The shape argument is "full" (default), "same" (the central part, the size of u) or "valid" (only the part computed without zero-padded edges). Inputs may be of different types, the output is double unless an input is single, and since R2026b sparse vectors are accepted. deconv divides polynomials, the inverse operation. conv2 is the 2-D version.
x = [1 2 3 4 5 6];
b = [1 1 1]/3;
disp(conv(x,b))
disp(conv(x,b,'same'))
disp(conv(x,b,'valid'))
[q,r] = deconv(conv([1 2 3],[1 1]),[1 1]); disp(q); disp(r)Octave 11.3.0 output:
0.3333 1.0000 2.0000 3.0000 4.0000 5.0000 3.6667 2.0000
1.0000 2.0000 3.0000 4.0000 5.0000 3.6667
2 3 4 5
1 2 3
0 0 0 0The first and last values of the full result are the edge effects of the zero padding; "valid" drops them.
Source: MATLAB documentation: conv
The convolution theorem
Convolution in time is multiplication in frequency: pad both vectors to the full length n = length(x)+length(b)-1, multiply the transforms and invert. Without padding you get circular convolution, which wraps around:
x = [1 2 3 4 5 6]; b = [1 1 1]/3;
n = length(x) + length(b) - 1;
z = real(ifft(fft(x,n).*fft(b,n))); % equals conv(x,b)
disp(z)
disp(real(ifft(fft([1 2 3 4]).*fft([1 0 0 1])))) % circularOctave 11.3.0 output:
0.3333 1.0000 2.0000 3.0000 4.0000 5.0000 3.6667 2.0000
3 5 7 5The first line equals conv(x,b). The second is a circular convolution of two length-4 vectors (1+4 = 5 is the wrapped first element, not the linear result 1). For long inputs the FFT route is faster than the direct sum; for short filters call conv.
Digital filtering with filter
y = filter(b,a,x)
y = filter(b,a,x,zi)
y = filter(b,a,x,zi,dim)
[y,zf] = filter(___)filter(b,a,x) applies the rational transfer function with numerator coefficients b and denominator coefficients a to x. If a(1) is not 1, the coefficients are normalised by it, so it must be nonzero. The output has the size of x; for a matrix it acts down the columns, and dim chooses another dimension (filter(b,a,x,[],2) filters rows). The difference equation is a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... - a(2)*y(n-1) - .... By default the delays start at zero, so the first outputs of a moving average are averages that include assumed zeros: the documentation’s 3-point average of [2 1 6 2 4 3] starts 0.6667 1.0000 3.0000. Pass zi (length max(length(a),length(b))-1) to continue from earlier data and take zf from a previous call to filter a long record in chunks.
A moving average is b = ones(1,w)/w with a = 1 (a finite impulse response, FIR); an exponential smoother needs a nonzero denominator (infinite impulse response, IIR):
x = [1 2 3 4 5 6];
disp(filter([1 1 1]/3,1,x))
disp(filter(0.2,[1 -0.8],[1 1 1 1 1]))Octave 11.3.0 output:
0.3333 1.0000 2.0000 3.0000 4.0000 5.0000
0.2000 0.3600 0.4880 0.5904 0.6723The moving average equals the first length(x) values of conv(x,b), as the first output shows; the second is the step response of the smoother y(n) = 0.2*x(n) + 0.8*y(n-1), which converges to 1. A causal filter delays the signal by about half the window; for zero delay, base MATLAB offers centred movmean(x,w) and smoothdata, while filtfilt is in the Signal Processing Toolbox.
For the heatflow signal, a 5-point moving average filter(ones(1,5)/5,1,x) attenuates the 3 Hz vibration (period of 3.3 samples) much more than the 1 Hz heater cycle (not run here). Check it by comparing fft amplitudes before and after.
Source: MATLAB documentation: filter
The boundary with Signal Processing Toolbox
| Base MATLAB | Signal Processing Toolbox |
|---|---|
|
Power spectral density ( |
|
Cross-correlation and spectral coherence, resampling ( |
|
Filter design ( |
The right column lists Signal Processing Toolbox functions, documented in the toolbox’s own reference (Signal Processing Toolbox); see Toolboxes and Add-On Products for how to check what is installed (ver, license).
Where the books are outdated
-
Spectrum functions. Older code that estimates spectra with
psdorspecgramshould move topwelchandspectrogram, their Signal Processing Toolbox successors. -
Corrupted listings. The Pérez López listings lost many
*and^signs, so formulas such asx.^2or2*pi*f*tcannot be trusted from that text. -
Smoothing. Moving averages are written today with
movmeanorsmoothdata; thefilterform remains right when you need a causal filter or chunked processing. -
Scaling.
ffthas no1/Nscaling andiffthas; amplitude plots need the division byNshown above.