fseriesexpo.m

非常好的数字处理教程

strGui.tDelay = ScaleIt(get(handles.uislTDelay,'Value'),0,4,8,0,T/2,T);                    % tDelay
% Fourier series form
t = get(handles.uirbFig1,'Value') + 2*get(handles.uirbFig2,'Value') + 4*get(handles.uirbFig3,'Value');
switch t
    case 1, strGui.fourierForm = 'Quadrature';
    case 2, strGui.fourierForm = 'Polar';
    case 4, strGui.fourierForm = 'Exponential';
    otherwise, error(sprintf('Unknown Fourier Form',t));
end
% Num Terms
sNTerms = get(handles.uitxNumTerms,'String');
sNTerms(1:10) = [];
strGui.N = str2num(sNTerms);
% Plot component ai, bi, etc.
t = get(handles.uirbPlotComponent1,'Value') + 2*get(handles.uirbPlotComponent2,'Value') + 4*get(handles.uirbPlotComponent3,'Value') + 8*get(handles.uirbPlotComponent4,'Value');
switch t
    case 1, strGui.plotComponent = 1;
    case 2, strGui.plotComponent = 2;
    case 4, strGui.plotComponent = 3;
    case 8, strGui.plotComponent = 4;
    otherwise, error(sprintf('Unknown component type %g',t))
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CaculateSeries %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function strData = CaculateSeries(strGui)

% get GUI parameters
waveType = strGui.waveType;
w = strGui.w;
T = 2*pi/w;
A = strGui.A;
sym = strGui.sym;
yOff = strGui.yOff;
tDelay = strGui.tDelay;
fourierForm = strGui.fourierForm;
N = strGui.N;
plotComponent = strGui.plotComponent;

%%%%%%%%%%%%%%%%%%%%%%%%%%% Calculate the time-domain vectors %%%%%%%%%%%%%%%%%%%%%
% create output vectors
t = linspace(0,2,400);
y0 = zeros(size(t));
y1 = y0;
% First create a triangle wave, regardless of what was requested
T1 = T*sym;  % Triangle changes direction at T1, which may be from 0 to T
m1 = -2/T1; b1 = 1; % line equations for first part
m2 = 2/(T-T1); b2 = 1-m2*T;
tt = mod(t-tDelay,T) + tDelay;
for i=1:length(t)
    if tDelay<=tt(i) && tt(i)<=T1+tDelay
        triangle(i) = m1*(tt(i)-tDelay) + b1;
    elseif T1+tDelay < tt(i) && tt(i) <= T+tDelay
        triangle(i) = m2*(tt(i)-tDelay) + b2;
    else
        error('Impossible time bracket')
    end
end
% Next change it into the requested type
switch strGui.waveType
    case 'Sin'
        for i=1:length(y0)
            theta = (1-triangle)*pi/2;
            y0 = A*cos(theta) +yOff;
        end
    case 'Square'
        if sym<=1e-6,
            y0 = ones(size(triangle));
        elseif sym>=1-1e-6
            y0 = -ones(size(triangle));
        else
            symOffset = (sym-1/2)*2;
            y0 = A*sign(triangle-symOffset) + yOff;
        end
        Ta = (m1*tDelay-b1+2*sym-1)/m1;
        Tb = (m2*tDelay-b2+2*sym-1)/m2;
        if w<=1e-6
            Ta=0; Tb=0;
        end
    case 'Triangle'
        y0 = A*triangle + yOff;
    otherwise
        error('Unknown waveType')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Calculate the c's %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create vectors
f = zeros(1,N+1);
c = f;
% determine c
switch strGui.waveType
    case 'Sin'
        tt1=m1*m2*pi*T*yOff;
        tt2=m2*cos(b1*pi/2)-m1*cos(pi/2*(b2+m2*T));
        tt3=-m2*cos(pi/2*(b1+m1*T1))+m1*cos(pi/2*(b2+m2*T1));
        tt4=m1*m2*pi*T;
        c(1)=(tt1+2*A*(tt2+tt3))/tt4;
        if w<=1e-5
            c(1)=y0(1);
        else
            for n=1:N
                tt1=-1/(n*T*w);
                tt2=pi/2;
                tt3=m1^2*pi^2-4*n^2*w^2;
                tt4=m2^2*pi^2-4*n^2*w^2;
                tt5=tt2*(b1+m1*T1);
                tt6=tt2*(b2+m2*T);
                tt7=tt2*(b2+m2*T1);
                tt8=j*m1*pi*cos(b1*tt2)-2*n*w*sin(b1*tt2);
                tt9=-j*m1*pi*cos(tt5)+exp(j*n*T1*w)*tt8+2*n*w*sin(tt5);
                tt10=-j*m2*pi*cos(tt6)+2*n*w*sin(tt6);
                tt11=j*m2*pi*tt3*cos(tt7)+tt4*tt9-2*n*w*tt3*sin(tt7);
                tt12=exp(-j*n*w*(T+tDelay))*yOff*(exp(j*n*T*w)-1);
                tt13=exp(j*n*w*(T1+2*tDelay))*tt3*tt10+exp(j*n*w*(T+2*tDelay))*tt11;
                tt14=2*A*n*w*tt13*exp(-j*n*w*(T+T1+3*tDelay));
                c(n+1)=tt1*j*(tt14/(tt3*tt4));
            end
        end
    case 'Triangle'
        if w<=1e-5
            c(1)=y0(1);
        else
            c(1)=A*(m2*T^2+2*b2*(T-T1)+2*b1*T1+(m1-m2)*T1^2)/(2*T)+yOff;
            for n=1:N
                tt1=j*n*w;
                tt2=exp(-tt1*(T+T1+tDelay))/(n^2*T*w^2);
                tt3=-exp(tt1*(T+T1))*(m1+tt1*b1);
                tt4=exp(tt1*T1)*(m2+tt1*(b2+m2*T));
                tt5=exp(tt1*T)*(m1-m2+tt1*(b1-b2+m1*T1-m2*T1));
                tt6=-tt1*exp(tt1*T1)*(exp(tt1*T)-1)*yOff;
                c(n+1)=tt2*(A*(tt3+tt4+tt5)+tt6);
            end
        end
    case 'Square'
        c(1)=A*(T+2*(Ta-Tb))/T+yOff; % this is really c0
        for n=1:N
            tt1=j*n*w;
            tt2=j*exp(-tt1*(T+2*Ta+Tb))/(n*T*w);
            tt3=-2*A*exp(tt1*(T+2*Ta));
            tt4=A*exp(tt1*(Ta+Tb));
            tt5=(A-yOff)*exp(tt1*(T+Ta+Tb));
            tt6=yOff*exp(tt1*(Ta+Tb));
            c(n+1)=tt2*(tt3+tt4+tt5+tt6);
        end
    otherwise
        error('unknown waveType')
end
% if elements approximately zero, make exactly zero
for n=1:N+1
    if abs(real(c(n)))<=1e-5
        c(n)=j*imag(c(n));
    end
    if abs(imag(c(n)))<=1e-5
        c(n)=real(c(n));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Synthesize f(t) %%%%%%%%%%%%%%%%%%%%%%%%%%%%
y1=y1+c(1);
for n=1:N
    y1=y1 + c(n+1)*exp(j*n*w*t) + conj(c(n+1))*exp(-j*n*w*t);
end
switch fourierForm
    case 'Quadrature'
        switch plotComponent
            case 1 %an
                yFreq=2*real(c);
            case 2 %bn
                yFreq=2*imag(c);
                yFreq(1)=0;
            otherwise
                error('bad plotComponent combination')
        end
    case 'Polar'
        switch plotComponent
            case 1 %An
                yFreq=2*abs(c);
            case 2 %theta
                yFreq=angle(c)*180/pi;
            otherwise
                error('bad plotComponent combination')
        end
    case 'Exponential'
        switch plotComponent
            case 1 %real
                yFreq=real(c);
            case 2 %imag
                yFreq=imag(c);
            case 3 %mag
                yFreq=abs(c);
            case 4 %angle
                yFreq=angle(c)*180/pi;
            otherwise
                error('bad plotComponent combination')
        end
end
if w<=1e-6
    xFreq=zeros(size(c));
else
    xFreq= linspace(0,N*w/(2*pi),N+1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%% Fill up output structure %%%%%%%%%%%%%%%%%%%%%%
strData.xTime = t;
strData.y0Time = y0;
strData.y1Time = y1;
strData.xFreq = xFreq;
strData.yFreq = yFreq;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ScaleIt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out=ScaleIt(sliderIn,sliderMin,sliderMean,sliderMax,outMin,outMean,outMax)
% ScaleIt provides a piecewise scaling of value SliderIn
% if SliderIn is between sliderMin and SliderMean it is scaled between outMin and outMean
% and similarly betwen Mean and Max
if sliderIn<sliderMin || sliderIn > sliderMax
    error('slider in outside min/max boundaries in ScaleIt')
end
if sliderIn < sliderMean
    tm = (outMean-outMin)/(sliderMean-sliderMin);
else
    tm = (outMax-outMean)/(sliderMax-sliderMean);
end
tb = outMean - (sliderMean * tm);
out = sliderIn * tm + tb;