Spatial filtering of cylindrical data

I have a set of dimensional measurements of a tire that are stored in a 2D array. Each column consists of 4096 measurements that are evenly spaced around the tire. The columns are also evenly spaced across the tire's surface, so the values in the array are the r coordinates, the rows are theta, and the columns are Z.

Assuming the radius of the tire is 1 foot, the circumference would be ~6.28 feet and the circumferential measurements would be 0.018 inch apart. The columns are about 1 mm apart (0.0394 inch), so the x and Y spatial resolutions are different.

I would like to high-pass and low-pass filter the data, but I am struggling to define an appropriate filter. I would probably filter the rows over a length of 6 inches (325 rows), but

1) how do I define such a filter in dimensional terms?

To remove any edge effects of the filtering, I can extend the row data by about 30% on both ends by adding data from the opposite ends of the array (which I will remove again after filtering).

Initially, I would like to only filter the data in the row direction, without a phase-shift, so filtfilt would seem appropriate to me.

2) Is that the right tool to use?

It would also be nice to be able to filter the entire surface (e.g., filter2), but I'm again not sure how to define such a 2D filter given the different spatial resolutions.

3) suggestions?

I do not have data available to pad the columns, so:

4) how do I deal with the edge effects of the filter on the rows?

I do have the signal processing toolbox and the image processing toolbox installed.

[Matt J]

the rows are theta, and the columns are Z.

Does that mean each row corresponds to a fixed theta, or does theta vary across each fixed row?

[Jim McIntyre]

Each row is a fixed theta: (rowNum - 1) * 2 * pi / 4096 radians.

Each column is a fixed lateral position: (columnNum - 1) mm

Answer 1

[William Rose]

@Jim McIntyre,

I think filtfilt() is a reasonable idea if you do what you said, which is add a wrap-around from the other part of the tire, at the beginning and end of the circle.

First we simulate some data.

R=304.8;                % tire radius (mm)
W=325;                  % tire width (mm)
N=4096;                 % number of angles
M=325;                  % number of z-values
theta=(0:N-1)*2*pi/N;   % angle (vector, radians)
c=R*theta;              % circumerential position (vector, mm)
z=linspace(-W/2,W/2,M); % axial position (vector, mm)

[C,Z]=meshgrid(z,theta); % size(C)=size(Z)=(4096,325)
Xz=1-(2*z/W).^2;         % tread height=parabolic function of z
% Xz=0 (fully worn) on the edges, =1 (not worn at all) in the middle.
Xt=1.5+cos(3*theta)+0.5*sin(360*theta);  % tread height as a function of theta
% Tread height has a large variation that repeats every 120 degrees around 
% the tire, and a smaller variation that repeats every 1 degree around. 
X=Xt'*Xz; % total tread height = product of circumferential and axial
figure
subplot(311), plot(z,Xz,'-r')
xlabel('Axial Position (mm)'); ylabel('Height Factor')
grid on; title("Tread Height vs Axial Pos.")
subplot(312), plot(c,Xt,'-r')
xlabel('Circumferential Position (mm)'); ylabel('Height Factor')
grid on; title("Tread Height vs Circum.")
subplot(313), surf(z,c,X,'EdgeColor','None');
xlabel('Axial'); ylabel('Circum (mm)'); zlabel('Tread')
colorbar; view(90,90); title('X=Tread Height')

(There appears to be aliasing, probably related to plot resolution, in the surface plot above, which causes some unexpected color changes in the horizontal direction. This aliasing does not happen when I run the script in regular Matab on my nortebook computer, rather than the Matlab Answers window.)

Make an augmented data array = [last half circumference; full tire circumference; first half circumference].

Xaug=[X(N/2+1:end,:);X;X(1:N/2,:)]; % size(Xaug)=8192x325

Filter the data going circumferentially around the tire with a zero-phase low-pass filter with a cutoff frequency corresponding to 150 mm (6").

dc =2*pi*R/N;               % circumferential sampling interval (mm)
fs=1/dc;                    % sampling rate (1/mm)
fco=1/150;                 % cutoff frequency (1/mm)
[b,a]=butter(2,fco/(fs/2)); % 2nd order lowpass Butterworth
Yaug=zeros(size(Xaug));     % allocate array
for j=1:M
    Yaug(:,j)=filtfilt(b,a,Xaug(:,j));
end
Y=Yaug(N/2+1:3*N/2,:);  % extract central part of filtered signal

Plot the filtered signal. Use "axis equal" to display with true aspect ratio.

figure
subplot(211)
surf(z,c,Y,'EdgeColor','None');
xlabel('Axial (mm)'); ylabel('Circumference (mm)');
title('Y=Circumferential Lowpass Filter'); 
axis equal; colorbar; view(90,90)

The lowpass filtered signal, Y, doesn't look much different than the plot of the original signal, X. (Aliasing makes them look diferent here, but the surface plots look quite simlar when I run it in regular Matab.) We can see the effect of lowpass filtering more clearly if we plot X versus circumference, on the axial centerline, for one third of the tire, before and after filtering.

subplot(212);
plot(c(1:1360),X(1:1360,(M+1)/2),'-r',c(1:1360),Y(1:1360,(M+1)/2),'-b')
legend('Raw','Filtered'); grid on; 
xlabel('Circumferential Position (mm)'); ylabel('Tread Height')
title('Tread Height along Axial Centerline for 1/3 Circumference')

[Jim McIntyre]

Thank you. This is extraordinarily helpful for my understanding of these kinds of filters. I expecially appreciate the effort you went to to document the parameters.

I extended what you wrote to also run a highpass filter and it leaves the "tread grooves" while removing the "non-uniformity."

The one question that remains is how to filter in two dimensions if I want to do that, while also preventing phase shift.

[William Rose]

@Jim McIntyre, youre welcome. Good job with the high pass filter.

For the 2D case, two options (among many) are

1. Use filtfilt() to apply a one-dimensional filter to each column of the 2D array, instead of applying it to each row. You could do this before or after filtering it in the circumferential direction as we did above.
2. Apply a 2D (spatial) filter to the data.

If you do (1), then you must think about the edges. The tread doesn't wrap around in the axial direction the way it does in the circumferential. (The tread part is not a torus, even though the tire is!) filtfilt() will handle the edges of a signal without you having to extend it at the edges, ands usually it give reasonable results. You can experiment with different ways of extending the edges to see which approach gives results that you deem acceptable. Then you clip off the etended edges of the filtered signal, of course.

If you do (2), you can use filter2() or conv2().

[William Rose]

@Jim McIntyre,

Here is a modified version of the earlier script. This version makes a different tread pattern, and shows how you can do 2D lowpass filtering with conv2(). As previously, I extend the tread surface circumferentially by wrapping around the tread from the other ends. Since I am filtering in 2D, I also extend the tread at the axial edges. I add a strip of points on each side, with vertical height=mean height of entire tread surface. The added strips have the same width as the moving average filter which I use to do the 2d lowpass filter. After filtering, I save the part of the extended, filtered surface that corresponds to the original tread area.

If the moving average filter dimensions are odd, then there is no phase shift. If the moving average filter dimensions are even, then there is a phase shift of 1/2 sample width. In this example, the moving average filter is square, but it doesn't have to be. A non-square filter would have different spatial cutoff frequencies in the circumferential and axial directions.

This script uses a flat moving average filter to do lowpass filtering with conv2(). You can also do other kinds of 2D lowpass or highpass filters. See imfilter(). And see here for discussion of conv2 versus filter2 versus imfilter.

[Jim McIntyre]

William,

Thanks! Your ideas have been really helpful. I took the code you gave me and modified it slightly below to give a slightly better representation of a tire. I then attempted to add both high and low pass filters, though I'm not certain that the 2D high pass filter is correct. I experimented with LC values from 1 to 150 and the results are not exactly what I was expecting.

Like you said earlier, the on-line graphics appear to show some aliasing, but when running it in "real" Matlab, the plots look better.

% Tire filtering, based on an answer from W. Rose 20250512
% See Matlab Answers post by @Jim McIntyre, 20250512.
% First we simulate some tire data
clear                                               % clear the workspace
Radius =     304.8;                                 % tire radius (mm)
Width =      325;                                   % tire width (mm)
CircSamps = 4096;                                   % number of angles
LatSamps =   325;                                   % number of z-values
Theta = (0 : CircSamps - 1) * 2 * pi / CircSamps;   % angle in circ direction (vector, radians)
Phi =   (0 : LatSamps - 1) * 2 * pi / LatSamps;     % angle in lat direction 
CircPos = Radius * Theta;                           % circumerential position (vector, mm)
LatPos = linspace(-Width/2, Width/2, LatSamps);     % axial position (vector, mm)
[C, Z] = meshgrid(LatPos, Theta);                   % size(C)=size(Z)=(4096,325)

% Xz = 0 (fully worn) on the edges, = 1 (not worn at all) in the middle.
% tread height is a parabolic function of z, with 5 tread ribs
Profile = 1 - (2 * LatPos / Width) .^ 2 + 0.5 * (cos(5 * Phi) < 0.7);  

% Tread height has a large variation that repeats every 120 degrees around 
% the tire, and a smaller variation that repeats every 1 degree around. 
CircVariation = 1.5 + cos(3 * Theta) + 0.5 * (sin(360 * Theta) < 0.8);  

% total tread height = product of circumferential and axial profiles
TreadRaw = CircVariation' * Profile + Radius; 
% Plot the simulated data
figure
subplot(411), plot(LatPos, Profile, '-r')
    xlabel('Axial Position (mm)'); ylabel('Height Factor')
    title("Tread Height vs Axial Pos.")
    grid on; 

subplot(412), plot(CircPos, CircVariation, '-r')
    xlabel('Circumferential Position (mm)'); ylabel('Height Factor')
    title("Tread Height vs Circum.")
    grid on; 

subplot(413), plot(CircPos, CircVariation, '-r')
    xlabel('Circumferential Position (mm)'); ylabel('Height Factor')
    title("Tread Height vs Circum. (detail)")
    xlim([0, 50]);
    grid on;

subplot(414), surf(LatPos, CircPos, TreadRaw, 'EdgeColor', 'None');
    xlabel('Axial'); ylabel('Circum (mm)'); zlabel('Tread')
    title('Raw Tread Height')
    colorbar; 
    view(90, 90); 

% (There appears to be aliasing, probably related to plot resolution, in the surface plot above, which causes some unexpected color changes in the horizontal direction.  This aliasing does not happen when I run the script in regular Matab on my nortebook computer, rather than the Matlab Answers window.)
% Make an augmented data array by wrapping data from both ends 
TreadAug = [TreadRaw(CircSamps / 2 + 1 : end, :); ...
            TreadRaw; 
            TreadRaw(1 : CircSamps / 2, :)]; % size(Xaug) = 8192 x 325
% Filter the data circumferentially around the tire 
dc = 2 * pi * Radius / CircSamps;   % circumferential sampling interval (mm)
fs = 1 / dc;                        % sampling rate (1/mm)
fco = 1 / 150;                      % cutoff frequency (1/mm)

% Build zero-phase low-pass and high-pass filters with cutoff frequency corresponding to 150 mm (6").
[bLow, aLow] =   butter(2, fco / (fs / 2), 'low');  % 2nd order lowpass Butterworth
[bHigh, aHigh] = butter(2, fco / (fs / 2), 'high'); % 2nd order highpass Butterworth

% allocate arrays
TreadFiltAugLow =  zeros(size(TreadAug));           
TreadFiltAugHigh = zeros(size(TreadAug));           

for j = 1 : LatSamps                % filter each column
    TreadFiltAugLow(:, j) =  filtfilt(bLow, aLow, TreadAug(:, j));
    TreadFiltAugHigh(:, j) = filtfilt(bHigh, aHigh, TreadAug(:, j));
end

% extract the central part of the filtered signals
TreadFiltLow =  TreadFiltAugLow(CircSamps / 2 + 1 : 3 * CircSamps / 2, :);   
TreadFiltHigh = TreadFiltAugHigh(CircSamps / 2 + 1 : 3 * CircSamps / 2, :); % + mean(TreadRaw,"all")
% Plot the filtered signal.  Use "axis equal" to display with true aspect ratio.
figure
subplot(311)
surf(LatPos, CircPos, TreadFiltLow, 'EdgeColor', 'None');
    xlabel('Axial (mm)'); ylabel('Circumference (mm)');
    title('Circumferential Lowpass Filter'); 
    axis equal; 
    colorbar; 
    view(90, 90)

subplot(312)
surf(LatPos, CircPos, TreadFiltHigh, 'EdgeColor', 'None');
    xlabel('Axial (mm)'); ylabel('Circumference (mm)');
    title('Circumferential Highpass Filter'); 
    axis equal; 
    colorbar; 
    view(90, 90)
% We can see the effect of lowpass filtering more clearly if we plot X versus circumference, on the axial centerline, for one third of the tire, before and after filtering.
subplot(313);
yyaxis left
plot(CircPos(1 : 1360), TreadRaw(1 : 1360, (LatSamps + 1) / 2),'-g', ...
     CircPos(1 : 1360), TreadFiltLow(1 : 1360, (LatSamps + 1) / 2), '-b')
    xlabel('Circumferential Position (mm)'); ylabel('Tread Height')
hold on
    yyaxis right
    plot(CircPos(1 : 1360), TreadFiltHigh(1 : 1360, (LatSamps + 1) / 2), '-r')
        ylim([-2,2]);
hold off
    legend('Raw','Low Pass','High Pass'); grid on; 
    xlabel('Circumferential Position (mm)'); ylabel('Tread Height')
    title('Tread Height along Axial Centerline for 1/3 Circumference')

% Filter the data in two directions
% Filter data: lowpass with square flat-top moving average.
dc = 2 * pi * Radius / CircSamps;   % circumferential sampling interval (mm)
dz = Width / LatSamps;              % axial sampling interval (mm)
Lc = 15;                            % filter size (mm)
nc = round(Lc / dc);                % circumferential extent of filter (samples)
nz = round(Lc / dz);                % axial extent of filter (samples)
BLow = ones(nc, nz) / (nc * nz);    % flat-top moving average
BHigh = -ones(nc, nz);
BHigh(ceil(end / 2), ceil(end / 2)) = -sum(BHigh,'all') - 1;
% BHigh(:, ceil(end / 2)) = (-sum(BHigh,'all') - 1) / nz;

% TreadAug2: pad shoulders by copying the overall mean to the left and on the right
borderStripe = mean(TreadRaw,"all") * ones(length(TreadAug), nz);
TreadAug2 = [borderStripe, TreadAug, borderStripe];

% Filter the augmented signal
TreadFiltAugLow = conv2(TreadAug2, BLow, "same");
TreadFiltAugHigh = conv2(TreadAug2, BHigh, "same");

% Extract the central part corresponding to original tread region
TreadFiltLow = TreadFiltAugLow(CircSamps / 2 + 1 : 3 * CircSamps / 2, nz + 1 : nz + LatSamps);
TreadFiltHigh = TreadFiltAugHigh(CircSamps / 2 + 1 : 3 * CircSamps / 2, nz + 1 : nz + LatSamps);

% Display the sizes of the different matrices, in the order that they were created
fprintf('Size: Tread=%dx%d, TreadAug=%dx%d, TreadAug2=%dx%d, TreadFiltAug2=%dx%d, Y=%dx%d.\n',...
    size(TreadRaw), size(TreadAug), size(TreadAug2), size(TreadFiltAugLow), size(TreadFiltLow))

Size: Tread=4096x325, TreadAug=8192x325, TreadAug2=8192x355, TreadFiltAug2=8192x355, Y=4096x325.

% Display the 2D filtered results
figure
subplot(311)
surf(LatPos,CircPos,TreadRaw,'EdgeColor','None');
xlabel('Axial (mm)'); ylabel('Circum (mm)'); zlabel('Tread')
colorbar; view(90,90); title('Raw Tread Height')
clim([305, 310]);           % set range for colorbar 
ylim([0 600]);              % 600 mm of circumference 

% Plot filtered signal
subplot(312)
surf(LatPos, CircPos, TreadFiltLow, 'EdgeColor', 'None');
xlabel('Axial (mm)'); ylabel('Circum (mm)')
title('2D Lowpass Filter'); 
clim([305, 310]);            
ylim([0, 600]);              
colorbar; view(90, 90)    

% Plot filtered signal
subplot(313)
surf(LatPos, CircPos, TreadFiltHigh, 'EdgeColor', 'None');
xlabel('Axial (mm)'); ylabel('Circum (mm)')
title('2D Highpass Filter'); 
clim([min(TreadFiltHigh,[],'all'), max(TreadFiltHigh,[],'all')]);            
ylim([0, 600]);              
colorbar; view(90, 90)    

[William Rose]

@Jim McIntyre, you're welcome. It looks to me like you're doing great with this.

Answer 2

[Matt J]

Using padarray you can pad circularly in theta and symmetrically in z and then filter with 'valid' shape flags. If you are filtering separably, you can also do fft-based convolution along the columns (in lieu of padding), which will be equivalent to a circulant filter in theta.

[Jim McIntyre]

Thank you. This is helpful in answering my question 4.