fs_fit.m

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

function [y,ty]=fs_fit(x,t,M,T)
% [y,ty]=fs_fit(x,t,N,T)
%
% Fourier-series fit. Fits a Fourier series to data
% samples not uniformly spaced in time.
% Inputs:
% x=vector of data samples.
% t=vector of corresponding time samples.
% M=number of frequencies in Fourier series, including the
%   dc term, the fundamental, and M-2 harmonics.
% T=time step (in units of t) of output y.
% Output:
% y=output vector with time step T. The length of y, and thus
%   the fundamental period of the Fourier series in units
%   of t, is ceil(max((t)-min(t))/T)*T. That is, the fit is
%   made over at least the duration of the input data. The
%   beginning of the period of y is at min(t).
% ty=vector of times corresponding with elements of y.
% Note: To test fs_fit using a regularly-sampled function,
% max(t)-min(t) should equal one complete period.

% Error checks etc.
if length(t)~=length(x)
    error('In fs_fit, t and x must have the same lenth');
elseif M<2
    error('In fs_fit, M must be at least 2.');
elseif max(t)-min(t)<3*T
    error('In fs_fit, duration of x is < 3T');
end
Nx=length(x);                       % Nx= # input samples
Ny=round((max(t)-min(t))/T)+1;      % Ny= length of output, y
G=zeros(Nx,2*M-1);                  % initialize G in (2.11)
G(:,1)=1/2;                         % first column of G
% The length of one period is Ny*T units of time, as stated above.
% Thus the arguments are in the form 2*pi*m*t(n)/(Ny*T). That is,
arg=2*pi*col_vec(t)/(Ny*T);         % vector of arguments; m=1
G(2:M,:)=cos(arg*[1:M-1]);          % columns 2-M in G
G(M+1:2*M-1,:)=sin(arg*[1:M-1]);    % columns M+1-2M-1 in G
c=(G'*G)\(row_vec(x)*G)';           % least-sq. coeff. (2.13)
c
ty=min(t)+[0:Ny-1]*T;               % time vector for y
arg2=2*pi*row_vec(ty)/(Ny*T);       % arguments for (2.23)
% The following expresses the Fourier series in (2.23).
y=c(1)+c(2:M)'*cos([1:M-1]'*arg2)+c(M+1:2*M-1)'*sin([1:M-1]'*arg2);