function [ dat ] = demod( sig, M, numchan, ts, sf, cf, fs  )
%[ dat ] = [ dat ] = demod( sig, M, numchan, ts, sf, cf, fs  
%   sig = received signal
%   M = bits per channel
%   N = Number of channels
%   ts = symbol duration
%   sf = safety factor (BW per channel = 1/(sf*ts)
%   cf = center frequency
%   fs = sample frequency



%First figure out the bandwidth for each signal

bwsig = sf/(ts);

chanbw = 2^M*bwsig;
totbw = numchan*chanbw;

str = sprintf('Total required bandwidth = %f', totbw);
disp(str)

T = 1/fs;

%Find channel center frequencies

%want ch1 + ch2 + ch3 + ch4 + ... + chn/n = cf
chanspace = totbw/numchan;

startbw = cf - totbw/2;

i = 1:numchan;
cbws = (i-1/2) * chanspace+startbw;

%Get FFT of incoming signal

symlen = ts/T;
siglen = length(sig)/symlen;

dat = -1 * ones(numchan, siglen);

numfreqs=2^M;
disp(symlen)
for i = 1:numchan
    stfreq = cbws(i) - numfreqs*bwsig/2;
    j = 1:numfreqs;
    cfreq = stfreq + (j-0.5)*bwsig;
    disp(cfreq)
    
    %Get energy in each of 2^M bins
    for j=0:siglen-1


        t = 0:T:ts-T;
        ecos = cos(2*pi*transpose(cfreq)*t);
        esin = cos(2*pi*transpose(cfreq)*t);

        symsig = sig(j*symlen+1:(j+1)*symlen);
        energy = zeros(2^M,1);
        for k = 1:2^M
           %energy(k) = sum( ( ecos(k,:).*symsig).^2 + (esin(k,:).*symsig).^2); 
           
           rp = sum(T*ecos(k,:).*symsig);
           ip = sum(T*esin(k,:).*symsig);
           energy(k) = rp.^2 + ip.^2;
           
        end

        [val, ind] = max(transpose(energy));
        dat(i, j+1) = ind-1;
    end
end
