(http://developer.mbed.org/handbook/Matlab-FIR-Filter)
320 samples at 48kHz for a sum of two sinusoids (1000Hz + 15000 Hz)
sample_rate = 48000;
nsamples = 256;
F = [1 15] * 1000;
A = [1 0.5];
% Time vector - use colon operator to generate integer vector of sample
% numbers
t = (0:nsamples-1) / sample_rate;
% Test signal - use matrix notation to compose it with single expression
signal = A * sin(2*pi*F'*t);
FIR filter
Assume a lowpass filter with cutoff frequency of 6 kHz.
The expectation is this should filter out the 15 kHz component from the test signal
% Choose filter cutoff frequency (6 kHz)
cutoff_hz = 6000;
% Normalize cutoff frequency (wrt Nyquist frequency)
nyq_freq = sample_rate / 2;
cutoff_norm = cutoff_hz / nyq_freq;
% FIR filter order (i.e. number of coefficients - 1)
order = 28;
% Create lowpass FIR filter through a direct approach: provide
% (normalized) cutoff frequency and filter order (assumed as known).
% fir1 takes care of designing the filter by imposing the constraints in
% the frequency domain and transforming back to time using a given window
% (the dafault used here is the Hamming window).
% For more advanced requirements see e.g. firpmord and firpm
% NOTE: fir1, firpmord and firpm all require Signal Processing Toolbox
fir_coeff = fir1(order, cutoff_norm);
% Analyse the filter using the Filter Visualization Tool
fvtool(fir_coeff, 'Fs', sample_rate)
% Filter the signal with the FIR filter
filtered_signal = filter(fir_coeff, 1, signal);
Plot the Signals
Also align filtered signal with original and discard transient samples (in this case the first order samples)
% Group delay as a scalar, in number of samples
group_delay = median(grpdelay(fir_coeff));
% Group delay in seconds
group_delay_s = group_delay / sample_rate;
% Plot the original signal...
figure(1)
plot(t, signal)
% ...and allow adding more plots
hold on
% Align and plot the filtered signal
% (On top of existing one)
plot(t-group_delay_s, filtered_signal, 'r-')
% Align and plot only the d6esired part of the filtered signal (discarding
% the transient)
plot(t(order:end)-group_delay_s, filtered_signal(order:end), ...
'g', 'LineWidth', 4);
grid on
hold off
