Spatial filtering of cylindrical data

Question

Jim McIntyre 2025-5-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2177096-spatial-filtering-of-cylindrical-data

评论： William Rose 2025-5-14

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.

2 个评论
显示无隐藏无

Matt J 2025-5-12

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 2025-5-12

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 2025-5-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2177096-spatial-filtering-of-cylindrical-data#answer_1565097

在 MATLAB Online 中打开

@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')

5 个评论
显示 3更早的评论隐藏 3更早的评论

William Rose 2025-5-13

tireTreadFilter2Db.m

@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 2025-5-14

在 MATLAB Online 中打开

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 2025-5-14

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

请先登录，再进行评论。

Answer 2

Matt J 2025-5-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2177096-spatial-filtering-of-cylindrical-data#answer_1565095

编辑：Matt J 2025-5-12

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jim McIntyre 2025-5-12

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

请先登录，再进行评论。

Spatial filtering of cylindrical data

2 个评论
显示无隐藏无

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Spatial filtering of cylindrical data

2 个评论 显示 无隐藏 无

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论