Differentiating in one direction using FFT2
显示 更早的评论
Starting from a univariateL-periodic function u, sampled at
evenly-spaced points in
, and stored in the vector u, the following function approximates the first derivative
, using FFT
function du = FFTDiff(u,n,L)
% Frequency vector
k = (2*pi/L) * [0:n/2-1, -n/2:-1]';
% Compute FFT
uHat = fft(u);
% Compute the derivative in the frequency domain, and back to physical domain
du = ifft( 1i * k .* uHat, 'symmetric');
end
First of all, if you see any alternative ways on how to write the function above, your feedback is welcome. I am now in the process of writing a function that does the following: starting from a bivariate, periodic function u, sampled at
evenly spaced points in
, I want to write a function that approximates
(hence the partial derivative with respect to x), using FFT2. The main problem I have is that I don't quite know how FFT2 outputs the wavenumbers
and
. It's a bit tricky to understand from the documentation, and I wonder if you have some ideas on how to do that.
3 个评论
Perhaps you have your reasons, but I thought I should point out that differentiating with FFTs is a highly inefficient thing to do. You should use diff() or gradient().
There might be cases where it makes sense if you are using some large and exotic differencing kernel, but that doesn't appear to be the case in the code you've shown.
Hi Daniele,
Can you show how this code works with a simple input? I'm running into an error.
Is the code supposed to work for L even and L odd?
L = 20; % period of function
l = 0:L; % n + 1 equally spaced points in [0,L]
n = numel(l) - 1; % n
u = sin(2*pi/20*l).'; % u
du = FFTDiff(u,n,L);
figure
plot(l,du,'-o')
function du = FFTDiff(u,n,L)
% Frequency vector
k = (2*pi/L) * [0:n/2-1, -n/2:-1]';
% Compute FFT
uHat = fft(u);
% Compute the derivative in the frequency domain, and back to physical domain
du = ifft( 1i * k .* uHat, 'symmetric');
end
Walter Roberson
2025-2-19
k = (2*pi/L) * [0:n/2-1, -n/2:-1]';
should probably be
k = (2*pi/L) * [0:n/2-1, -n/2:-1/2]';
采纳的回答
更多回答(1 个)
How about this? It generalizes your original function to let you differentiate along any specified dimension dim for any nD array u -
function du = FFTDiff(u,L,dim)
arguments
u double {mustBeNonempty}
L (1,1) double
dim (1,1) double = find(size(u)>1,1) %default to first non-singleton dimension
end
n=size(u,dim);
% Frequency vector
k = (2*pi/L) * ifftshift( (0:n-1)-ceil((n-1)/2) );
e=ones(1,ndims(u)); e(dim)=n;
k=reshape(k,e);
% Compute FFT
uHat = fft(u,[],dim);
% Compute the derivative in the frequency domain, and back to physical domain
du = ifft( 1i * k .* uHat, [],dim,'symmetric');
end
2 个评论
Walter Roberson
2025-2-18
I think your code might possibly have some trouble of u is empty. In that case, size(u, find(size(u)>1,1)) becomes size(u,[]) which returns [] . But [] cannot be stored in something size (1,1)
If you are successful in storing the [] into dim, then n=size(u,dim) would return empty, and I suspect that would lead to problems.
Catalytic
2025-2-18
Yes, probably. I fixed it.
类别
在 帮助中心 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!