amod.m

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

function [y, t] = amod(x, Fc, Fs, method, opt1, opt2, opt3, opt4, opt5)
%AMOD   Analog signal modulation.
%       Y = amod(X, Fc, Fs, METHOD, OPT1, OPT2, OPT3) modulates the
%       message signal X with carrier frequency Fc (in Hz) and sample
%       frequency Fs (Hz) using the modulation scheme in string variable
%       METHOD. The time interval between two successive points in X 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. For best result, Fs > Fc is required (by
%       Nyquist theory).
%
%       METHOD is a string, which can be one of the following:
%       'amdsb-tc'  Amplitude modulation, double sideband with transmission
%                   carrier
%       'amdsb-sc'  Amplitude modulation, double sideband suppressed carrier
%       'amssb'     Amplitude modulation, single side band suppressed carrier
%       'qam'       Quadrature amplitude modulation
%       'fm'        Frequency modulation
%       'pm'        Phase modulation
%
%       Use AMOD(METHOD) to view the help for a specific method.
%
%       See also: ADEMOD, AMODCE, DMOD, DMODCE.

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

opt_pos = 5;
if nargin < 1
    disp('Usage: Y=AMOD(X, Fc, Fs, METHOD, OPT1, OPT2, OPT3)');
    return;
elseif isstr(x) & nargin == 1
    % help lines for individual modulation method.
    method = deblank(x);
    if strcmp(method, 'am')
        method = 'amdsb-tc';
    end;
    addition = 'See also COMDEMOD, DMOD, DDEMOD, AMODCE, ADEMODCE.';
    callhelp('amod.hlp',method, addition);
    return;
elseif nargin < 3
    disp('Usage: Y=AMOD(X, Fc, Fs, METHOD, OPT1, OPT2, OPT3) for analog modulation');
    return;
elseif nargin < 4    
    method = 'am';
end

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

method = lower(method);

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

%begin the calculation.
if findstr(method, 'amdsb-sc')
    t = (0:1/Fs:((size(x,1)-1)/Fs))';
    t = t(:, ones(1, size(x, 2)));
    y = x .* cos(2 * pi * Fc * t + ini_phase);
elseif findstr(method, 'amdsb-tc') | strcmp(method, 'am')
    if nargin < opt_pos
        opt1 = -min(min(x));
    end;
    t = (0:1/Fs:((size(x, 1)-1)/Fs))';
    t = t(:, ones(1, size(x, 2)));
    y = (x + opt1) .* cos(2 * pi * Fc * t + ini_phase);
elseif findstr(method, 'amssb')
    t = (0:1/Fs:((size(x,1)-1)/Fs))';
    t = t(:, ones(1, size(x, 2)));
    if nargin < opt_pos
        % do the Hilbert transform in frequence domain.
        if findstr(method, '/up')
            y = x .* cos(2 * pi * Fc * t + ini_phase) - ...
                imag(hilbert(x)) .* sin(2 * pi * Fc * t + ini_phase);    
        else
            y = x .* cos(2 * pi * Fc * t + ini_phase) + ...
                imag(hilbert(x)) .* sin(2 * pi * Fc * t + ini_phase);    
        end;
        return;
    elseif nargin <= opt_pos
        [num, den] = hilbiir(1/Fs);
    else
        num = opt1;
        den = opt2;
    end;
    y = x;
    for  i = 1 : size(x, 2)
        y(:,i) = filter(num, den, y(:,i));
    end;
    if findstr(method, '/up')
        y = x .* cos(2 * pi * Fc * t + ini_phase) - ...
            y .* sin(2 * pi * Fc * t + ini_phase);
    else
        y = x .* cos(2 * pi * Fc * t + ini_phase) + ...
            y .* sin(2 * pi * Fc * t + ini_phase);
    end; 
elseif strcmp(method, 'qam')
    if floor(c/2) ~= c/2
        error('QAM in amod takes even column matrix only.')
    end;
    [len_x, wid_x] = size(x);
    if ceil(wid_x/2) ~= wid_x/2
        error('Input signal for AMOD QAM must have even number of columns.')
    end;
    t = (0:1/Fs:((size(x,1)-1)/Fs))';
    t = t(:, ones(1, wid_x/2));
    y = x(:, 1:2:wid_x) .* cos(2 * pi * Fc * t + ini_phase) + ...
        x(:, 2:2:wid_x) .* sin(2 * pi * Fc * t + ini_phase);
elseif strcmp(method, 'fm')
    if nargin >= opt_pos
        x = opt1 * x;
    end;
    t = (0:1/Fs:((size(x,1)-1)/Fs))';
    t = t(:, ones(1, size(x, 2)));
    x = 2 / Fs * pi * x;
    x = [zeros(1, size(x, 2)); cumsum(x(1:size(x,1)-1, :))];
    y = cos(2 * pi * Fc * t + x + ini_phase);
elseif strcmp(method, 'pm')
    if nargin >= 6
        x = opt1 * x;
    end;
    t = (0:1/Fs:((size(x, 1)-1)/Fs))';
    t = t(:, ones(1, size(x, 2)));
    y = cos(2 * pi * Fc * t + x + ini_phase);
else
    % digital modulation, format changed to be
    % Y = AMOD(X, Fc, Fs, METHOD, Fd, M, OPT1, OPT2, OPT3);
    if nargin < 5
        M = max(max(x)) + 1;
        M = 2^(ceil(log(M)/log(2)));
        M = max(2, M);
    else
        M = opt1;
    end;

    if nargin < 5
        ini_phase = 0;
    else
        ini_phase = opt1;
    end;

    if strcmp(method, 'ask')
        x = (x - (M-1)/2)*2/(M-1);
        [y, t] = amod(x, Fc, Fs, 'amdsb-sc', ini_phase);
    elseif strcmp(method, 'psk')
        x = 2 * pi * x / M;
        [y, t] = amod(x, Fc, Fs, 'pm', ini_phase);
    elseif findstr(method, 'qask')
        if strcmp(method, 'qask')
            %plain square QASK style.
            [s1, s2] = qaskenco(M);
        elseif findstr(method, 'qask/cir')
            if nargin < opt_pos+1
                disp('Usage: Y = AMOD(X, Fc, Fs, ''qask/cir'', NIC, AIC, PIC).')
                return
            elseif nargin < opt_pos+2
                opt3 = 0;
            end;
            NIC = M;
            AIC = ini_phase;
            PIC = opt2
            s1 = apkconst(NIC, AIC, PIC);
            s2=imag(s1);
            s1 = real(s1);
            ini_phase = opt3;
            M = sum(NIC);
        else
            s1 = M;
            s2 = ini_phase;
            if nargin < opt_pos
                disp('Usage: Y = AMOD(X, Fc, Fs, ''qask/arb'', In_Phase, Quad).')
            elseif nargin < opt_pos + 1
                opt2 = 0;
            end
            ini_phase = opt2;
            M = length(s1);
        end;
        if min(min(x)) < 0
            error('The input signal for MODULPE has illegal element.');
        elseif max(max(x)) > M
            error('The input signal for MODULPE has element larger than given limit.');
        end;
        x = x(:) + 1;
        y = amod([s1(x) s2(x)], Fc, Fs, 'qam', ini_phase);
    elseif strcmp(method, 'fsk')
        % M has been calculated
        % Fd = ini_phase (opt1)
        if nargin < opt_pos
            Fd = Fs / M;
        else
            Fd = opt1;
        end;
        % Tone = opt2
        if nargin < opt_pos + 1
            Tone = 2 * Fd / M;
        else
            Tone = opt2;
        end;
        % ini_phase = opt3
        if nargin < opt_pos + 2
            ini_phase = 0;
        else
            ini_phase = opt3;
        end;
        % modulation map part.
        FsDFd = Fs/Fd;
        offset = 0;
        if length(Fd) > 1
           offset = rem(rem(Fd(2), FsDFd) + FsDFd, FsDFd);
        end;
        %To make the modulation counting only the very first number of each digit
        %transfer, unquote the following five lines.
        %x = x(:)';
        %y = x(1+offset:FsDFd:length(x));
        %y = y(ones(1, FsDFd), :);
        %y = y(:);
        %x(1+offset:length(x)) = y(1:length(x)-offset)';
        x = x(:) * Tone;
        % analog modulation.
        y = [];
        len_x = length(x);
        mod_b = 1;
        mod_e = FsDFd;
        if offset > 0
            mod_e = offset;
        end;
        while mod_b <= len_x
            y = [y; amod(x(mod_b:min(len_x, mod_e)), Fc, Fs, 'fm', ini_phase)];
            mod_b = mod_b + FsDFd;
            mod_e = mod_e + FsDFd;
        end;
        t = (0:1/Fs:((size(x,1)-1)/Fs))';
        t = t(:, ones(1, size(x, 2)));
    elseif strcmp(method, 'msk')
        %msk is a special case of fsk, with M=2; Tone = Fd.
        M = 2;
        if nargin < opt_pos - 1
            Fd = Fs/2;
        else
            Fd = ini_phase;
        end;
        if nargin < opt_pos
            ini_phase = 0;
        else
            ini_phase = opt1;
        end;
        % it is a special case of fsk with M = 2 and Tone = Fd. 
        y = amod(X, Fc, Fs, 'fsk', ini_phase, 2, Fd, Fd, ini_phase);
    else
        %not machtch any one of them. Display error message.
        disp(['Modualtion method ', method, ' is not a legal option in function AMOD.'])
        disp('The current available mehtods are:');
        disp('  amdsb-tc Amplitude modulation, double sideband with transmission carrier.')
        disp('  amdsb-sc Amplitude modulation, double sideband supressed carrier.')
        disp('  amssb Amplitude modulation, single side band supressed carrier.')
        disp('  qam Quadrature Amplitude modulation.')
        disp('  pm Phase modulation.')
        disp('  ask Amplitude shift keying modulation.')
        disp('  psk Phase shift keying modulation.')
        disp('  qask Quadrature amplitude shift-keying modulation, square constellation.')
        disp('  qask/cir Quadrature amplitude shift-keying modulation, circle constellation.')
        disp('  qask/arb Quadrature amplitude shift-keying modulation, user defined constellation.')
        disp('  fsk Frequency shift keying modulation.')
        disp('  msk Minimum shift keying modulation.')
    end;
end;
t = t(:, 1);
if (r==1)
    y = y.';
    t = t.';
end;

% end amod.m