Mismatch in Banded Toeplitz Matrix Method for FIR Filter Output

Hi MATLAB Community,

I am working on a script that designs and applies a simple FIR low-pass filter. The goal is to compare the filter output using three approaches:

1. The filter function (reference method).
2. A banded Toeplitz matrix of the signal. (BTsig)
3. A banded Toeplitz matrix of the filter coefficients. (BTfilt)

Please find the full code at :

https://github.com/cyber-g/Matlab-Filter-Playground/blob/main/filter_experiment.m

While the results from the filter function and the BTfilt coefficients match closely, the output from the BTsig method does not align with the other two. I know there is an issue on line 54:

Tx_tmp = conj(tril(toeplitz(x(:,nrand)))); % Toeplitz matrix

but I cannot find the right processing to do.

Thank you in advance for your help!

Answer 1

[AH]

I think that this Tx_tmp = tril(toeplitz(x(:,nrand),[x(1,nrand),zeros(1,length(b)-1)])); % Toeplitz matrix will solve the problem.

% filter_experiment.m
% This script designs a simple FIR low-pass filter. 
% It generates a complex Gaussian test signal, applies the filter, and compares 
% the output using two methods: 
% 1) A banded Toeplitz matrix of the signal.
% 2) A banded Toeplitz matrix of the filter coefficients.
% The script visualizes the results and computes the Mean Squared Error (MSE) 
% for both methods to evaluate their accuracy.
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also: filter
% Author: Germain Pham
% email: [e-mail]
% April 2025; Last revision: 
% 
% Special thanks to Denis Gilbert, who provided the original code for "M-file Header Template"
% https://fr.mathworks.com/matlabcentral/fileexchange/4908-m-file-header-template
% which is adapted here for the purpose of this example.
% 
% Special thanks to the unkown author of:
% https://www.dsprelated.com/freebooks/filters/Matrix_Filter_Representations.html
% 
%------------- BEGIN CODE --------------

% Design a simple FIR low pass filter
% Define the filter specifications
Fs = 1000; % Sampling frequency
Fc = 100;  % Cut-off frequency
N  = 50;   % Filter order

% Generate the low pass FIR filter coefficients using Hamming window
b = fir1(N, Fc/(Fs/2), hamming(N+1));
b = b(:); % Ensure b is a column vector

% Generate a test signal
Nsamples = 2000; % Number of samples
Nrandomize = 10; % Number of realizations
x = randn(Nsamples, Nrandomize) + 1j * randn(Nsamples, Nrandomize); % Complex Gaussian noise
x = x ./ (ones(Nsamples,1)*sqrt(var(x))); % Normalize the signal

% Apply the filter to the test signal
y_ref_filter = filter(b, 1, x); % Filter the signal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Method 1: Using the banded Toeplitz SIGNAL matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the Tril-toeplitz matrix of the signal
Tx = zeros(Nsamples,length(b),Nrandomize);
for nrand = 1:Nrandomize
    Tx_tmp = tril(toeplitz(x(:,nrand),[x(1,nrand),zeros(1,length(b)-1)])); % Toeplitz matrix
    Tx(:,:,nrand) = Tx_tmp(:,1:length(b));     % Keep only the first N columns corresponding to the filter length
end

% Compute the filter output using matrix multiplication
y_matrix_x = squeeze(pagemtimes(Tx,repmat(b(:),[1 1 Nrandomize]))); % Matrix multiplication

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Method 2: Using the banded Toeplitz FILTER matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Tb = tril(toeplitz([b;zeros(length(x)-length(b),1)])); % Toeplitz matrix for filter coefficients

% Compute the filter output using matrix multiplication
y_matrix_b = Tb * x; % Matrix multiplication


% Plots
nexttile
plot(real([y_ref_filter(:,1) y_matrix_x(:,1) y_matrix_b(:,1)]))
nexttile
plot(imag([y_ref_filter(:,1) y_matrix_x(:,1) y_matrix_b(:,1)]))
legend('Reference Filter', 'Matrix Method (x)', 'Matrix Method (b)')
title('Filter Output Comparison')
xlabel('Sample Index')
ylabel('Amplitude')
grid on

% Errors 
MSE_x = sqrt(mean(abs(y_ref_filter - y_matrix_x).^2)); % Mean Squared Error for x
MSE_b = sqrt(mean(abs(y_ref_filter - y_matrix_b).^2)); % Mean Squared Error for b

disp(['MSE for banded Toeplitz SIGNAL: ', num2str(MSE_x)]);

MSE for banded Toeplitz SIGNAL: 1.3166e-16  1.3074e-16  1.2771e-16  1.3133e-16  1.2544e-16  1.2786e-16  1.2381e-16   1.308e-16  1.3608e-16  1.3029e-16

disp(['MSE for banded Toeplitz FILTER: ', num2str(MSE_b)]);

MSE for banded Toeplitz FILTER: 2.5775e-17   2.507e-17  2.3583e-17  2.3436e-17  2.8353e-17  2.4714e-17  2.6961e-17  2.1099e-17  2.8609e-17  2.2841e-17