% This function generate the sample sequences from the signal and the
% specified MWC system

function [Samples, TimeAxis, A, Sig, Nyq_TimeAxis] = MWC_SubNyquistSample(Signal,MWC,handles)

% The procedure of generating the lowrate samples, in numerical simulation,
% consists of the following stages:
% 1. generate the input signal at the Nyquist rate
% 2. generate the sequences z[n], by appropriate mixing and filtering.
% These sequences are decimated to rate f_s.
% 3. Let q=floor(Tp/Ts); A=S*F*D ; if q>1 then update A according to the
% expansion factor q
% 4. Mix the sequences z[n] according to the matrix A, and decimate the
% results to 1/Ts rate

% 1: generate the signal at the Nyquist rate (or higher)

% We would like to avoid noninteger decimations, therefore we generate the
% signal on a time-axis that is at least at the Nyquist resolution,
% and that allows for integer decimations in the sequel

fmax = getfield(Signal.Structure, 'fmax');
fnyq = 2*fmax;
M = getfield(MWC, 'M');%, ceil(fnyq/fp))
m = getfield(MWC, 'm');%, ceil(fnyq/fp))
fp = fnyq/M;
Ts = getfield(MWC, 'Ts');
fs = 1/Ts;
decfactor = ceil(fnyq/fp); % this is the decimation ratio from the Nyquist-rate to fp
fgrid = decfactor*fp;  % fgrid>=fnyq, and both fgrid/fp and fgrid/fs are integers

% get the required number of samples
NumSamples = str2double(get(handles.editNumSamples,'String'));
Npts = NumSamples*decfactor*2;  % we take a factor of 2.. just to be safe

%Set to Empty
CleanBar(handles);
%Show Bar
ShowBar(handles);
set(handles.uneditStatus,'String','Status: Busy, Generating signal');
CleanBar(handles)
[SignalAmp, NoiseAmp, Nyq_TimeAxis] = GenerateSignalValues(Signal,1/fgrid,Npts,handles);
Sig = SignalAmp + NoiseAmp;
% 2: Generate z[n]
L = ceil(fnyq/fp);
if (mod(L,2)==0)
    L=L+1;
end   % L is the number of specturm slices, namely the length of z[n]
L0 = (L-1)/2;
%z = zeros(L,length(Nyq_TimeAxis));

% prepare a polyphase filter H(f)
% Analog filter (polyphase for speed)
dec = fgrid/fp;
h_analog = fir1(5000, 1/dec);
lenh = length(h_analog);
polylen = ceil(lenh / dec);
h_pad = [h_analog zeros(1,polylen*dec-lenh)];
polyh = reshape(h_pad,dec,polylen);
CleanBar(handles);
ShowBar(handles);
set(handles.uneditStatus,'String','Status: Busy, Generating samples');

%% 3: Prepare sensing matrix
S = getfield(MWC, 'SignPatterns');
DFTM = fft(eye(M));
F = [ conj(fliplr(DFTM(:,2:(L0+1))))    DFTM(:,1:(L0+1)) ];
theta = exp(-j*2*pi/L);
np = 1:L0;
nn = (-L0):1:-1;
dn = [   (1-theta.^nn)./(j*2*pi*nn)      1/L    (1-theta.^np)./(j*2*pi*np)  ];
D = diag(dn);
A = conj(S*F*D);

z=0;
for Zind = 1:L
    slice = Zind - L0 - 1;
    % x(t) -> X(f)
    % X(f-D) is shifting the spectrum to the right by D Hz
    % X(f-D) -> x(t)exp(j 2 pi D t)
    % for slice=-L0, we want to shift the spectrum by L0*fs to the
    % right... this explains the minus
    Sig_shifted = Sig.*exp(-j*2*pi*slice*fp*Nyq_TimeAxis);
    [SamplesZ, Fp_TimeAxis] = PolyFilterDecimate(Sig_shifted,Nyq_TimeAxis,polyh,lenh);
    if (Zind==1)
        y_tilde=zeros(m,length(SamplesZ));
    end
    % 4: Sensing and interpolation to 1/Ts rate
    %y_tilde = C*z;
    y_tilde = y_tilde + A(:,Zind)*SamplesZ;
    updateProgBar(Zind,handles,L*2);
end

TimeAxis = Fp_TimeAxis;
Samples = y_tilde;
