integral.m

英文书《Digital Signal Processing with Examples in MATLAB》附带的MATLAB实例

function y=integral(x,fs)
% y=integral(x,fs)
% 
% y =frequency-domain integration of x
%    in a column vector or array of column vectors.
% fs =sampling frequency (samples/s).
% Unlike other methods (trapezoidal, splines, etc.), y is
% exact, provided x has no frequencies at or above fs/2.
% Note: if x is an array, each column of x is integrated.
% Note: if length(col. of x) is odd, the last element is
%       ignored and length(col. of y) =length(col. of x)-1.

T=1/fs;                             %T =time step (s)
[N,nc]=size(x);
if N==1,                            % if x is a row vector,
    x=x';                           %make x a column vector
    [N,nc]=size(x);                 %N=length of each vector
end
% Check for errors.
if N<4,
   error('vector(x) must have at least 4 elements.');
elseif rem(N,2)==1,
    N=N-1;                          %ignore last sample if N is odd
end
y=zeros(N,nc);                      %output vector or array
k=[0:N-1]';                         %k =sample #
m=[1:N/2-1]';                       %m =frequency indices
for col=1:nc,                       %integrate each col. of x
    X=fft(x(1:N,col));              %X =DFT{column vector}
    Im=imag(X(m+1));                %see slide on integration
    V0=-N*T*sum(Im./m)/pi;
    Vm=-N*T*j*X(m+1)./(2*pi*m);
    V=[V0; Vm; 0; rev(conj(Vm))];
    c1=k*T*X(1)/N;                  %dc component
    y(:,col)=c1+ifft(V);
end