【正文】
(39。Value for dataout_type must be a 1, 2 or 339。)。end%Determine Number of points (next power of 2), frequency increment%and Nyquist frequencyN = 2^nextpow2(max(size(datain)))。df = 1/(N*dt)。Nyq = 1/(2*dt)。%Save frequency arrayiomega_array = 1i*2*pi*(Nyq : df : Nyqdf)。iomega_exp = dataout_type datain_type。%Pad datain array with zeros (if needed)size1 = size(datain,1)。size2 = size(datain,2)。if (Nsize1 ~= 0 amp。amp。 Nsize2 ~= 0)if size1 size2datain = vertcat(datain,zeros(Nsize1,1))。elsedatain = horzcat(datain,zeros(1,Nsize2))。endend%Transform datain into frequency domain via FFT and shift output (A)%so that zerofrequency amplitude is in the middle of the array%(instead of the beginning)A = fft(datain)。A = fftshift(A)。%Convert datain of type datain_type to type dataout_typefor j = 1 : Nif iomega_array(j) ~= 0A(j) = A(j) * (iomega_array(j) ^ iomega_exp)。elseA(j) = plex(,)。endend%Shift new frequencyamplitude array back to MATLAB format and%transform back into the time domain via the inverse FFT.A = ifftshift(A)。datain = ifft(A)。%Remove zeros that were added to datain in order to pad to next%biggerst power of 2 and return dataout.if size1 size2dataout = real(datain(1:size1,size2))。elsedataout = real(datain(size1,1:size2))。endreturn