ademod.m

数字通信第四版原书的例程

function x  = ademod(y, Fc, Fs, method, opt1, opt2, opt3, opt4, opt5)
%ADEMOD Demodulates a modulated signal.
%       X = ADEMOD(Y, Fc, Fs, METHOD, OPT1, OPT2, OPT3, OPT4, OPT5)
%       demodulates a modulated signal Y with carrier frequency Fc (Hz) and
%       sample frequency Fs (Hz) using the demodulation method provided in
%       METHOD. The time interval between two successive points in Y is 1/Fs.
%       When Fs is a two element vector, the second element specifies the
%       initial phase in the carrier signal modulation. The unit of the
%       initial phase is rad. NUM and DEN are the numerator and
%       denominator of the low-pass filter respectively. The sample time of
%       the filter is 1/Fs. When NUM is zero, empty, or not given, this
%       function uses the default low-pass filter which is
%       [NUM, DEN] = BUTTER(5, Fc*2/Fs).
%
%       Use ADEMOD(METHOD) to view the help for a specific method.
%
%       METHOD is a string, which can be one of the following:
%       'amdsb-tc'  Amplitude demodulation, double sideband with transmission
%                   carrier.
%       'amdsb-sc'  Amplitude demodulation, double sideband suppressed
%                   carrier.
%       'amssb'     Amplitude demodulation, single side band suppressed
%                   carrier.
%       'qam'       Quadrature amplitude demodulation.
%       'fm'        Frequency demodulation.
%       'pm'        Phase demodulation.
%
%       See also: AMOD.
%

%       Wes Wang 9/29/95
%       Copyright (c) 1995-96 by The MathWorks, Inc.
%       $Revision: 1.1 $  $Date: 1996/04/01 17:51:30 $
%       This function calls ademod.help for further help.

opt_pos = 5;
if nargin < 1
    disp('Usage: X = ADEMOD(Y, Fc, Fs, METHOD, OPT1, OPT2, OPT3, OPT4, OPT5)');
    return;
elseif isstr(y) & nargin == 1
    % help lines for individual modulation method.
    method = deblank(y);
    if strcmp(method, 'am')
        method = 'amdsb-tc';
    end;
    addition = 'See also: ADEMOD, DMOD, DDEMOD, AMODCE, ADEMODCE.';
    callhelp('ademod.hlp',method, addition);
    return;
end;
if nargin < 4
    disp('Usage: X = ADEMOD(Y, Fc, Fs, METHOD, INI_PHASE, NUM, DEN, OPT1, OPT2, OPT3');
    return;
end;

method = lower(method);

[r, c] = size(y);
if r*c == 0
    y = [];
    return;
end;
if (r == 1)
    y = y(:);
    len = c;
else
    len = r;
end;
pi2 = 2 * pi;

if length(Fs) < 2
    ini_phase = 0;
else
    ini_phase = Fs(2);
    Fs = Fs(1);
end;

%begin the processing
%process amdsb-sc and amdsb-tc, am is the same as amdsb-sc.
if ~isempty(findstr(method, 'amdsb')) | ...
        ((length(method)==2) & strcmp(method(1:2), 'am')) | ...
        ~isempty(findstr(method, 'am/'))
    if findstr(method, 'amdsb-sc')
        if nargin < opt_pos+1
            [num,den] = butter(5, Fc*2/Fs);
        else
            num = opt1;
            den = opt2;
        end;
    else
        if nargin < opt_pos + 2
            [num,den] = butter(5, Fc*2/Fs);
        else
            num = opt2;
            den = opt3;
        end;
    end;
    if findstr(method, 'costas');
        % costas method. reference to SIMULINK block diagram for the algorithm.

        %pre-process the filter.
        if abs(den(1)) < eps
            error('First denominator filter coeficient must be non-zero.');
        else
            num = num/den(1);
            if (length(den) > 1)
                den = - den(2:length(den)) / den(1);
            else
                den = 0;
            end;
            num = num(:)';
            den = den(:)';
        end;
        len_den = length(den);
        len_num = length(num);

        y = y * 2;
        for ii = 1:size(y,2)
            z1 = zeros(length(den), 1);
            z2 = z1;
            z3 = z1;
            s1 = zeros(len_num, 1);
            s2 = s1;
            s3=s1;
            intgl = 0;

            for i = 1:size(y,1);
                %begining from the integration value.
                fm_out_i = cos(pi2 * intgl + ini_phase);
                fm_out_q = -sin(pi2 * intgl + ini_phase);
                s1 = [y(i, ii)*fm_out_i; s1(1:len_num-1)];
                s2 = [y(i, ii)*fm_out_q; s2(1:len_num-1)];
                x(i, ii) = num * s1 + den * z1;
                z1 = [x(i, ii); z1(1:len_den-1)];
                z2 = [num*s2+den*z2; z2(1:len_den-1)];
                s3 = [z1(1)*z2(1); s3(1:len_num-1)];
                z3 = [num*s3+den*z3; z3(1:len_den-1)];
                intgl = rem((z3(1) + Fc)/Fs + intgl, 1);
            end;
        end;
    else
        t = (0 : 1/Fs :(len-1)/Fs)';
        t = t(:, ones(1, size(y, 2)));
        x = y .* cos(pi2 * Fc * t + ini_phase);
        for i = 1 : size(y, 2)
            x(:, i) = filter(num, den, x(:, i)) * 2;
        end;
    end;
    if findstr(method, 'amdsb-tc')
        if nargin >= opt_pos
            x = x - opt1;
        else
            x = x - mean(x);
        end;
    end;
elseif findstr(method, 'amssb')
    if nargin < opt_pos+1
        [num, den] = butter(5, Fc * 2 / Fs);
    else
        num = opt1;
        den = opt2;
    end;
    t = (0 : 1/Fs :(len-1)/Fs)';
    t = t(:, ones(1, size(y, 2)));
    x = y .* cos(pi2 * Fc * t + ini_phase);
    for i = 1 : size(y, 2)
        x(:, i) = filter(num, den, x(:, i)) * 2;
    end;
elseif findstr(method, 'qam')
    if nargin < opt_pos+1
        [num, den] = butter(5, Fc * 2 / Fs);
    else
        num = opt1;
        den = opt2;
    end;
    t = (0 : 1/Fs :(len-1)/Fs)';
    t = t(:, ones(1, size(y, 2)));
    x = [y .* cos(pi2 * Fc * t + ini_phase) y.*sin(pi2*Fc*t + ini_phase)];
    wid_x = size(x, 2);
    tmp = [1:wid_x/2;wid_x/2+1:wid_x];
    tmp = tmp(:);
    x = x(:,tmp);
    for i = 1 : wid_x
        x(:, i) = filter(num, den, x(:, i)) * 2;
    end;
elseif ~isempty(findstr(method, 'fm')) | ~isempty(findstr(method, 'pm'))
    is_fm = findstr(method, 'fm');
    if nargin < opt_pos+1
        [num, den] = butter(5, Fc * 2 / Fs);
    else
        num = opt1;
        den = opt2;
    end;
    if isempty(num)
        [num, den] = butter(5, Fc * 2 / Fs);
    end;
    if nargin > opt_pos+1
        sen = opt3;
    else
        sen = 1;
    end;
    ini_phase = ini_phase + pi/2;
    %pre-process the filter.
    if abs(den(1)) < eps
        error('First denominator filter coeficient must be non-zero.');
    else
        num = num/den(1);
        if (length(den) > 1)
            den = - den(2:length(den)) / den(1);
        else
            den = 0;
        end;
        num = num(:)';
        den = den(:)';
    end;
    len_den = length(den);
    len_num = length(num);

    x = y;
    y = 2 * y;
    for ii = 1 : size(y, 2)
        z1 = zeros(length(den), 1);
        s1 = zeros(len_num, 1);
        intgl = 0;

        memo = 0;
        for i = 1:size(y, 1)
            %start with the zero-initial condition integer.
            vco_out = cos(pi2 * intgl + ini_phase);
            s1 = [y(i, ii) * vco_out; s1(1:len_num-1)];
            tmp = num * s1 + den * z1;
            z1 = [tmp; z1(1:len_den-1)];
            intgl = rem(((tmp*sen + Fc)/ Fs + intgl), 1);
            if is_fm
                x(i, ii) = tmp;
            else
                memo = memo + tmp * sen / Fs;
                x(i, ii) = memo * pi2;
            end;
        end;
    end;
    x = x;
else
    disp(['Method ', method, ' is not a legal option in function ADEMOD.']);
end;