function [x,x_soft] = ddemod_self(y, Fc, Fd, Fs, method, M, opt1, opt2, opt3, opt4)
%DDEMOD
%
%WARNING: This is an obsolete function and may be removed in the future.
% Please use MODEM.PAMDEMOD, MODEM.QAMDEMOD, MODEM.GENQAMDEMOD,
% FSKDEMOD, MODEM.PSKDEMOD, or MODEM.MSKDEMOD instead.
% Z = DDEMOD(Y, Fc, Fd, Fs, METHOD...) demodulates the digital modulated
% signal Y with carrier frequency Fc (Hz), sample frequency Fd (Hz), and
% computation sampling frequency Fs (Hz). The ratio Fs/Fd must be a
% positive integer. For information about METHOD and subsequent
% parameters, and about using a specific demodulation technique,
% type one of these commands at the MATLAB prompt:
%
% FOR DETAILS, TYPE DEMODULATION TECHNIQUE
% ddemod ask % M-ary amplitude shift keying demodulation
% ddemod psk % M-ary phase shift keying demodulation
% ddemod qask % M-ary quadrature amplitude shift keying
% % demodulation
% ddemod fsk % M-ary frequency shift keying demodulation
% ddemod msk % Minimum shift keying demodulation
% ddemod sample % Downsample Y without demodulating
%
% For baseband simulation, use DDEMODCE. To plot signal constellations,
% use MODMAP.
%
% See also DMOD, AMOD, ADEMOD, DMODCE, DDEMODCE, DEMODMAP, MODMAP,
% EYEDIAGRAM, SCATTERPLOT.
% Copyright 1996-2007 The MathWorks, Inc.
% $Revision: 1.1.6.5 $
% warnobsolete(mfilename, 'Please use MODEM.PAMDEMOD, MODEM.QAMDEMOD, MODEM.GENQAMDEMOD, FSKDEMOD, MODEM.PSKDEMOD, or MODEM.MSKDEMOD instead.');
swqaskenco = warning('off', 'comm:obsolete:qaskenco');
swapkconst = warning('off', 'comm:obsolete:apkconst');
swademod = warning('off', 'comm:obsolete:ademod');
swdemodmap = warning('off', 'comm:obsolete:demodmap');
opt_pos = 7; % position of 1st optional parameter
if nargin < 1
feval('help', 'ddemod');
return;
elseif isstr(y)
method = lower(deblank(y));
if length(method) < 3
error('Invalid method option for ddemod.');
end
if nargin == 1
addition = ['See also DMOD, AMOD, ADEMOD, DMODCE, DDEMODCE, DEMODMAP, MODMAP,',...
'\r EYEDIAGRAM, SCATTERPLOT.'];
addition = sprintf(addition);
if method(1:3) == 'qas'
callhelp('ddemod.hlp', method(1:4), addition);
else
callhelp('ddemod.hlp', method(1:3), addition);
end
else
warning('Wrong number of input variables. Use MODMAP to plot constellations.');
end
return;
end
if nargin < 4
disp('Usage: Z=DDEMOD(Y, Fc, Fd, Fs, METHOD, M, OPT1, OPT2, OPT3, OPT4) for passband demodulation');
return;
elseif nargin < opt_pos - 2
if nargout < 1
method = 'eye';
else
method = 'sample';
end
end
method = lower(method); % findstr is case sensitive
if length(Fs) > 1
ini_phase = Fs(2);
Fs = Fs(1);
else
ini_phase = 0; % default initial phase
end
if length(Fd) > 1
offset = Fd(2);
Fd = Fd(1);
else
offset = 0; % default timing offset
end
if ~isfinite(Fs) | ~isreal(Fs) | Fs<=0
error('Fs must be a positive number.');
elseif ~isfinite(Fd) | ~isreal(Fd) | Fd<=0
error('Fd must be a positive number.');
else
FsDFd = Fs/Fd; % oversampling rate
if ceil(FsDFd) ~= FsDFd
error('Fs/Fd must be a positive integer.');
end
end
if ~isreal(offset) | ceil(offset)~=offset | offset<0 | offset>=FsDFd
error('OFFSET must be an integer in the range [0, Fs/Fd).');
end
if length(Fc) ~= 1 | ~isfinite(Fc) | ~isreal(Fc) | Fc <= 0
error('Fc must be a positive number. For baseband demodulation, use DDEMODCE.');
elseif Fs/Fc < 2
warning('Fs/Fc must be much larger than 2 for accurate simulation.');
end
if (nargin >= opt_pos & isempty(findstr(method, '/arb')) & ...
isempty(findstr(method, '/cir')) & ...
(length(M) ~= 1 | ~isfinite(M) | ~isreal(M) | M <= 0 | ceil(M) ~= M))
error('Alphabet size M must be a positive integer.');
end
if isempty(y)
x = [];
return;
end
[r, c] = size(y);
if r == 1
y = y(:);
len_y = c;
else
len_y = r;
end
if rem(len_y, FsDFd) ~= 0
error('Number of samples in y must be an integer multiple of Fs/Fd.');
elseif ~isreal(y)
error('Input Y must be real.');
end
if strncmpi(method, 'ask', 3)
% M has to be specified in 'ask' option.
if nargin < opt_pos - 1
error('Not enough input variables for DDEMOD.');
end
index = ([0 : M - 1] - (M - 1) / 2) * 2 / (M - 1);
if nargin < opt_pos + 1
% In digital demodulation, integrator replaced LPF.
num = 1;
den = 1;
else
num = opt1;
den = opt2;
end
if findstr(method, '/cos')
y = ademod(y, Fc, [Fs, ini_phase], 'amdsb-sc/costas', num, den);
else
y = ademod(y, Fc, [Fs, ini_phase], 'amdsb-sc', num, den);
end
% Integrate to remove double freq component and replicate average
% over symbol
sizey = size(y);
y = integ(y, FsDFd, offset);
y = repmat(y(:), 1, FsDFd);
y = reshape(y.', sizey(1), sizey(2));
if findstr(method, '/eye')
ddemod(y, Fc, [Fd, offset], [Fs, ini_phase], 'eye');
end
if findstr(method, '/sca')
ddemod(y, Fc, [Fd, offset], [Fs, ini_phase], 'sca');
end
if findstr(method, '/nomap')
x = y;
else
x = demodmap(y, [Fd, offset], Fs, 'ask', M);
end
elseif strncmpi(method, 'fsk', 3)
if nargin < opt_pos
Tone = Fd;
else
Tone = opt1;
end
if findstr(method,'/nomap')
warning(sprintf(['The option ''/nomap'' does not apply to FSK demodulation.\n',...
' The function will proceed ignoring the ''/nomap'' switch.']));
end
%calculate the correlation of fsk.
[len_y, wid_y] = size(y);
z = [-(M-1):2:(M-1)] * Tone * pi / Fs;
z = [ones(len_y, 1)]*z;
z = cumsum(z);
t = [0 : 1/Fs : 1/Fd-1/Fs]';
t = t(:, ones(1, M));
symbol_period=1/Fd;
%leave space for x
x = y([offset+1 : FsDFd : len_y], :);
[len_x, wid_x] = size(x);
if findstr(method, '/eye')
t1 = [0 : FsDFd-1]/Fs;
t1 = t1 + offset/Fs;
clf;
plot([min(t1), max(t1), max(t1)], [-1/2, NaN, 1])
axis([min(t1) max(t1), -1/2, 1]);
hold on
end
for i = 1 : wid_x
comp_low = 1;
if offset <= 0
comp_upp = FsDFd;
else
comp_upp = offset;
end
for k = 1 : len_x
if findstr(method, '/nonc')
z_temp = cos((t+(k-1)*symbol_period)*2*pi*Fc + z(1:FsDFd,:));
zz_temp = sin((t+(k-1)*symbol_period)*2*pi*Fc + z(1:FsDFd,:));
else
% column i of z_temp represents correlation sequence for ith symbol
z_temp = cos((t + (k-1)*symbol_period)*2*pi*Fc + z(comp_low:comp_upp,:) + ini_phase);
end
if findstr(method, '/eye')
if findstr(method, '/nonc')
cor1 = y(comp_low:comp_upp, i)*ones(1,M) .* z_temp(1:comp_upp-comp_low + 1, :);
cor2 = y(comp_low:comp_upp, i)*ones(1,M) .* zz_temp(1:comp_upp-comp_low + 1, :);
corr = cumsum(cor1).^2 + cumsum(cor2).^2;
else
corr = y(comp_low:comp_upp, i)*ones(1,M) .* z_temp(1:comp_upp-comp_low+1,:);
corr = cumsum(corr);
end
corr = corr / max(max(max(abs(corr))), eps);
if k == 1
plot(t1(FsDFd-size(corr, 1)+1 : FsDFd)', corr);
else
plot(t1(1:size(corr, 1))', corr)
end
corr = corr(size(corr, 1), :);
else
if findstr(method, '/nonc')
cor1 = y(comp_low:comp_upp, i)*ones(1,M) .* z_temp(1:comp_upp-comp_low + 1, :);
cor2 = y(comp_low:comp_upp, i)*ones(1,M) .* zz_temp(1:comp_upp-comp_low + 1, :);
corr = sum(cor1).^2 + sum(cor2).^2;
else
corr = y(comp_low:comp_upp, i)*ones(1,M) .* z_temp(1:comp_upp-comp_low+1,:);
corr = sum(corr);
end
end
[tmp_v, tmp_i] = max(corr);
x(k, i) = tmp_i - 1;
% ------------ 计算软判决结果 -------------
% 大于1,硬判决为1,小于1硬判决为0
x_soft(k,i) = corr(2)/corr(1) ;
% ------------ 计算软判决结果 -----------