function [y,freq,c]=exzfft_ma(x,fe,fs,nfft,D)

nt=length(x);           % ¼ÆËã¶ÁÈëÊý¾Ý³¤¶È
fi=fe-fs/D/2;           % ¼ÆËãϸ»¯½ØֹƵÂÊÏÂÏÞ
fa=fi+fs/D;             % ¼ÆËãϸ»¯½ØֹƵÂÊÉÏÏÞ
na=round(0.5 * nt/D+1); % È·¶¨µÍͨÂ˲¨Æ÷½ØֹƵÂʶÔÓ¦µÄÆ×ÏßÌõÊý
% ƵÒÆ
n=0: nt-1;              % ÐòÁÐË÷ÒýºÅ
b=n*pi* (fi+fa)/fs;     % ÉèÖõ¥Î»ÐýתÒò×Ó
y=x.*exp(-1i*b);        % ½øÐÐƵÒÆ 
b= fft(y, nt);          % FFT
% µÍͨÂ˲¨ºÍϲÉÑù
a(1: na) =b(1: na);     % È¡ÕýƵÂʲ¿·ÖµÄµÍƵ³É·Ö
a(nt-na+2 : nt) =b(nt-na+2 : nt); % È¡¸ºÆµÂʲ¿·ÖµÄµÍƵ³É·Ö
b=ifft(a, nt);          % IFFT 
c=b(1 : D: nt);         % ϲÉÑù
% Çóϸ»¯ÆµÆ×
y=fft(c, nfft) * 2/nfft;% ÔÙÒ»´ÎFFT
y=fftshift(y);          % ÖØÐÂÅÅÁÐ
freq=fi+(0:nfft-1)*fs/D/nfft; % ƵÂÊÉèÖÃ

figure
plot(freq,20*log10(y))
grid on