Can MATLAB do a 2-D FFT on a uniform parallelogram grid of measured points?

67 次查看（过去 30 天）
Jerry Guern
Jerry Guern 2024-11-18，23:16
编辑： Paul about 7 hours 前
I have some measurements in slant coordinates, and I need to project them into the ground plane. When I just project the slant grid into ground coords, it ends up a uniform parallelogram grid. In order to form an image, I need to Fourier Transform my parallelogram grid of measured values into ground freq domain, the Inv FT that into ground time domain. The code below does this successfully, using a toy function f(x,y) = cos(3x-4y) in place of measurements. But in this code, I calculated non-fast Four Transform coefficients one at a time.
EDIT: I want to emphasize, this code does what I need, calculating the varialbe 'pfd' in Figure 41, I just need to know if some MATLAB function can somehow do it with a Fast Fourier Transform.
dim1 = 29; % Resolution of image we're creating
dim2 = 31;
% This block just creates xs & ys, coords of a random parallelogram of
% measured sample points.
theta1 = rand * pi/4;
theta2 = theta1 * (.5 + rand);
M = [2*(1+rand)/3 * [cos(theta1) sin(theta1)]; 2*(1+rand)/3 * [-sin(theta2) cos(theta2)]];
[Xs, Ys] = meshgrid(linspace(-3,3,300),linspace(-3,3,300));
gpt = M * [Xs(:)';Ys(:)']; %These points sprawl out all around the image area.
% Find the measured points that fall within the image.
keep = gpt(1,:)>-1/(2*dim1) & gpt(1,:)<1-1/(2*dim1) & gpt(2,:)>-1/(2*dim2) & gpt(2,:)<1-1/(2*dim2);
xs = gpt(1,keep);
ys = gpt(2,keep);
figure(1); scatter(xs,ys,.5); title('measured points within image')
% Create image pixels
[Xq, Yq] = meshgrid(linspace(0,1-1/(dim1),dim1),linspace(0,1-1/(dim2),dim2));
%Ratio of measured points within image to image pixels.
oversample = sum(keep) / (dim1 * dim2)
% Wave numbers for Fourier Transform
kx = [0 [1:floor(dim1/2)] [-floor(dim1/2) : -1]];
ky = [0 [1:floor(dim2/2)] [-floor(dim2/2) : -1]]';
% For this toy problem, use a closed form function in place of measured
% data. Calclate this function first for our image points.
td = cos(3*Xq-4*Yq);
% Calculate the frequency domain, from our td above and for our toy
% measured data.
jfd = zeros(size(td));
nufd = zeros(size(td));
for j = 1:dim1
for k = 1:dim2
% Calculate the non-fast FFT, just for comparison
jfd(k,j) = sum(sum(td .* exp(-1i*2*pi*(kx(j)*Xq+ky(k)*Yq))));
% Calculate same terms, but from our measured samples.
pfd(k,j) = sum(cos(3*xs-4*ys) .* exp(-1i*2*pi*(kx(j)*xs+ky(k)*ys)));
end
end
pfd = pfd / oversample;
figure(39); surf(Xq,Yq,td); title('orig function');
figure(40); surf(Xq,Yq,ifft2(jfd)); title('ifft2 of non-fast FT for comparison');
figure(41); surf(Xq,Yq,real(ifft2(pfd))); title('real of ifft2 of pfd');
So, is there any MATLAB function that do the same thing, calculate that pfd, but using an FFT?
I thought at one time that nufftn() did this, but I honestly couldn’t make any sense of the vague documentation or the comments in the code.
  2 个评论
Paul
Paul 2024-11-19，3:22
Hi Jerry,
If I'm following correctly, it seems that xs and ys define 1956 pairs of non-gridded points (xs and ys are both 1 x 1956). Is that the intent of the code? If so, my reading of nufftn is the "time" variables don't have to be equally spaced in each dimension, but they still have to be gridded, i.e., they can't be scattered points. In other words, we can have
t1 = sort(rand(1,N));
t2 = sort(rand(1,M));
and the function to be evaluated would be something like
z = cos(t1) + sin(t2');
Same thing for the frequency variables, i.e., they define a non-uniformly spaced grid (in general).
I'd try running some code here, but Answers doesn't seem to executing code right now (at least not for me).
Jerry Guern
Jerry Guern about 19 hours 前
Thanks for the reply. xs and ys will always have the same length, but that length will be different in each run, because they're on a randomly oriented grid . But it is always a uniform grid, albeit a slanted parallelogram grid. If you run the code, Figure 1 will show you what I mean. I tried to use nufftn() for this, but I couldn't make any sense of the documentation or code comments.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Matt J
Matt J about 17 hours 前
编辑：Matt J about 17 hours 前
I need to Fourier Transform my parallelogram grid of measured values into ground freq domain, the Inv FT that into ground time domain.
I don't know what "ground time/freq" domain means.
I think if you just pretend the samples are Cartesian and take the FFT and IFFT, you will obtain the non-Fourier domain result sampled on your original parallelgoram. You can then just use griddata or scatteredInterpolant to reample that, if needed.
  2 个评论
Jerry Guern
Jerry Guern about 16 hours 前
I'm creating an image in ground plane coordinates. When I said ground time/freq domain, the time domain would be the values of the ground plane image I'm trying to create, and the freq domain is the fft of that. I said ground to contrast it with the slant plane, the coordionates I'm projecting into ground plane. I don't know what you meant by what you said. Could you show me in code?

请先登录，再进行评论。

Chris Turnes
Chris Turnes about 16 hours 前
Yes, nufftn can compute those two loops:
jfd2 = reshape(nufftn(td, { Yq(:,1), Xq(1,:) }, { ky kx }), dim2, dim1);
pfd2 = reshape(nufftn(cos(3*xs-4*ys), [ys(:) xs(:)], { ky kx }), dim2, dim1) / oversample;
If you compute norm(jfd - jfd2) and norm(pfd - pfd2) you'll find they match to within around 1e-12, 1e-13 in value.
  1 个评论
Paul
Paul about 7 hours 前
编辑：Paul about 7 hours 前
@Chris Turnes
If I'm understanding pfd2 correctly ....
The 1D vectors xs and ys define pairs of points that are not gridded in x-y space. The sample points to nufftn uses [ ] to indicated that the sample points are given in a list. The query points input uses {} to indicate to nufftn that the query points are gridded. Is that correct?
Is the reshape really needed for pfd2?
Do you think some additional examples at nufftn would be helpful?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by