nufft在matlab中源代码

以下是MATLAB中nufft函数的源代码：

function [y] = nufft(x, om, N, K, n_shift, iflag)
%NUFFT  NUFFT algorithm (using convolution via FFT)
%  [y] = nufft(x, om, N, K, n_shift, iflag)
%  Input
%    x:      input signal (# of samples x # of channels)
%    om:     frequency locations (# of samples x 2)
%               k = om(:,1) * 2*pi / N; l = om(:,2);
%    N:      signal length
%    K:      # of Fourier modes computed (K>=N)
%    n_shift:  0 or 1 (default)
%              perform fftshift or not in the output
%    iflag:  sign of the imaginary unit in the exponential
%              if iflag>0,  e^(iwt) = cos(wt) + i*sin(wt)
%              if iflag<0,  e^(iwt) = cos(wt) - i*sin(wt)
%  Output
%    y:      output spectrum (# of modes x # of channels)
%
%  If you use this code, please cite the following paper:
%  Y. C. Eldar, A. J. Feuer, and G. D. Forney, Jr.,
%  "Optimal tight frames and quantum measurement," IEEE Trans. Inform. Theory, vol. 48, no. 3, pp. 599-610, Mar. 2002.
%
%  The NUFFT algorithm is based on the following paper:
%  A. Dutt and V. Rokhlin, "Fast Fourier transforms for nonequispaced data," SIAM J. Sci. Comput., vol. 14, no. 6, pp. 1368-1393, Nov. 1993.
%
%  The FFTSHIFT function is used to match the output of MATLAB's FFT function.
%
%  Written by Alex Pothen (alex.pothen@gmail.com), Sep 2013.
%  Modified by Yonina Eldar (yonina.eldar@gmail.com), Oct 2013.

if nargin < 6, iflag = 1; end
if nargin < 5, n_shift = 1; end

if iflag > 0
    sign = 1;
else
    sign = -1;
end

M = size(om,1);

% Compute Fourier transform of convolution kernel
kk = (-K/2):(K/2-1);
vk = exp(sign*2*pi*i*kk/N).';
NN = (-N/2):(N/2-1);
n_shift = N*mod(n_shift,2); % make sure n_shift is either 0 or N
h = zeros(N,1);
h(NN + n_shift + 1) = sqrt(N) * vk .* (abs(NN + n_shift) < K/2);
h = ifftshift(h);
h = fft(h);

% Compute convolution of signal with kernel
y = zeros(K,size(x,2));
for ii=1:size(x,2)
    y(:,ii) = h .* fftshift(fft(ifftshift(x(:,ii))));
end

% Compute Fourier transform at om frequencies
y = y(mod(om(:,2)),:);
y = y(om(:,1)+1,:);

% Perform fftshift if required
if n_shift
    y = fftshift(y,1);
end

这是一个比较基础的NUFFT实现，将输入信号x与卷积核h卷积，然后在频域对结果进行采样以得到输出y。