%
% OFDM signal generation.
%
%   Takes a text file as input and generates an audio file that is
%   modulated with N-QAM modulation and transmitted as a series of OFDM
%   packets.  Each packet contains a synchronization signal, channel
%   estimation signal, and the encoded data.
%
%   Usage
%      [ret] = pulsegen(in_file, type, out_file, [fs, sync, fft_size, 
%         sym_per_packet, cp_length])
%
%   Parameters
%      ret - OPTIONAL - writes the function return value or error code
%          if the function does not complete successfully.
%          Errors are:
%             0 - SUCCESS
%             1 - INVALID_MOD_TYPE
%             2 - INVALID_ARGS
%             3 - INPUT_FILE_ERROR
%      in_file - string that specifies a ASCII text file to encode
%      type - the constellation type.  Options are:
%          BPSK, QPSK, 8QAM, 16QAM, 32QAM, 64QAM
%      out_file - output file that the encoded data is written to.
%                 The file name is concatenated with '.wav'
%      fs - OPTIONAL - audio file samping rate in Hz. Default is 8 kHz.
%      sync - OPTIONAL - a column vector that containts a time-domain 
%         sync code.  This defaults to a 32 length sine wave code.
%      fft_size - OPTIONAL - the size of the FFT.  Default is 64.
%      sym_per_packet - OPTIONAL - the number of symbols per OFDM packet.
%         The default is 14.  If you increase this number you have less
%         overhead from syncing and channel estimation but you risk your
%         channel response changing so that the channel estimation no
%         longer corrects the symbols properly.
%      cp_length - OPTIONAL - length of the cyclic prefix.  Default is 8.
%         A longer prefix allows more leeway in packet sync but adds some
%         overhead in the time domain waveform.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change Log
% ----------------
% Caleb Kemere Feb., 2001 (created original version)
% Modified, Thomas Bruns March, 2001
% Modified, Jiwon Seo, Feb. 2008
% Modified, Andrew Price, Feb. 2009
%    - turned into function
%    - added support for generating multiple modulation types
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [varargout] = pulsegen(in_file, type, out_file, varargin)

   % check to make sure the required arguments are there
   if (nargout == 1)
       DISPLAY = 0;
   elseif (nargout == 0)
       DISPLAY = 1;
   else
       disp('ERROR: Incorrect number of output arguments.');
   end
   if (nargin < 3)
       if(DISPLAY)
          disp('ERROR: Not enough arguments supplied.');
       else
          varargout{1} = 2;
       end
       return;
   end
   
   % figure out what type of modulation was requested and initialize the
   % carrier location lookup table.
   switch (type)
       case 'BPSK'
           bitsPerCarrier = 1;
           carrierLookup = zeros(2^bitsPerCarrier, 2);
           carrierLookup(1,:) = [1 0]; %0
           carrierLookup(2,:) = [-1 0]; %1
       case 'QPSK'
           bitsPerCarrier = 2;
           carrierLookup = zeros(2^bitsPerCarrier, 2);
           carrierLookup(1,:) = [1 1]; %00
           carrierLookup(2,:) = [-1 1]; %01
           carrierLookup(3,:) = [1 -1]; %10
           carrierLookup(4,:) = [-1 -1]; %11
       case '8QAM'
           bitsPerCarrier = 3;
           carrierLookup = zeros(2^bitsPerCarrier, 2);
           carrierLookup(1,:) = [1 1]; %000
           carrierLookup(2,:) = [0 1]; %001
           carrierLookup(3,:) = [1 0]; %010
           carrierLookup(4,:) = [-1 1]; %011
           carrierLookup(5,:) = [0 -1]; %100
           carrierLookup(6,:) = [-1 -1]; %101
           carrierLookup(7,:) = [1 -1]; %110
           carrierLookup(8,:) = [-1 0]; %111
       case '16QAM'
           bitsPerCarrier = 4;
           carrierLookup = zeros(2^bitsPerCarrier, 2);
           carrierLookup(1,:) = [1/3 1/3]; %0000
           carrierLookup(2,:) = [1/3 1]; %0001
           carrierLookup(3,:) = [1 1/3]; %0010
           carrierLookup(4,:) = [1 1]; %0011
           carrierLookup(5,:) = [1/3 -1/3]; %0100
           carrierLookup(6,:) = [1/3 -1]; %0101
           carrierLookup(7,:) = [1 -1/3]; %0110
           carrierLookup(8,:) = [1 -1]; %0111
           carrierLookup(9,:) = [-1/3 1/3]; %1000
           carrierLookup(10,:) = [-1/3 1]; %1001
           carrierLookup(11,:) = [-1 1/3]; %1010
           carrierLookup(12,:) = [-1 1]; %1011
           carrierLookup(13,:) = [-1/3 -1/3]; %1100
           carrierLookup(14,:) = [-1/3 -1]; %1101
           carrierLookup(15,:) = [-1 -1/3]; %1110
           carrierLookup(16,:) = [-1 -1]; %1111
       case '32QAM'
           bitsPerCarrier = 5;
           carrierLookup = zeros(2^bitsPerCarrier, 2);
           carrierLookup(1,:) = [-1 1/4]; %00000
           carrierLookup(2,:) = [0 -3/4]; %00001
           carrierLookup(3,:) = [0 1/4]; %00010
           carrierLookup(4,:) = [1 1/4]; %00011
           carrierLookup(5,:) = [1 -1/4]; %00100
           carrierLookup(6,:) = [0 3/4]; %00101
           carrierLookup(7,:) = [0 -1/4]; %00110
           carrierLookup(8,:) = [-1 -1/4]; %00111
           carrierLookup(9,:) = [-1/2 3/4]; %01000
           carrierLookup(10,:) = [-1/2 -1/4]; %01001
           carrierLookup(11,:) = [1/2 3/4]; %01010
           carrierLookup(12,:) = [1/2 -1/4]; %01011
           carrierLookup(13,:) = [1/2 -3/4]; %01100
           carrierLookup(14,:) = [1/2 1/4]; %01101
           carrierLookup(15,:) = [-1/2 -3/4]; %01110
           carrierLookup(16,:) = [-1/2 1/4]; %01111
           carrierLookup(17,:) = [-3/4 -1/2]; %10000
           carrierLookup(18,:) = [1/4 -1/2]; %10001
           carrierLookup(19,:) = [-3/4 1/2]; %10010
           carrierLookup(20,:) = [1/4 1/2]; %10011
           carrierLookup(21,:) = [3/4 1/2]; %10100
           carrierLookup(22,:) = [-1/4 1/2]; %10101
           carrierLookup(23,:) = [3/4 -1/2]; %10110
           carrierLookup(24,:) = [-1/4 -1/2]; %10111
           carrierLookup(25,:) = [1/4 1]; %11000
           carrierLookup(26,:) = [-3/4 0]; %11001
           carrierLookup(27,:) = [1/4 0]; %11010
           carrierLookup(28,:) = [1/4 -1]; %11011
           carrierLookup(29,:) = [-1/4 -1]; %11100
           carrierLookup(30,:) = [3/4 0]; %11101
           carrierLookup(31,:) = [-1/4 0]; %11110
           carrierLookup(32,:) = [-1/4 1]; %11111
       case '64QAM'
           bitsPerCarrier = 6;
           carrierLookup = zeros(2^bitsPerCarrier, 2);
           % -- TODO --
       otherwise
           if (DISPLAY)
               disp('ERROR: Unrecognized modulation type.');
           else
               varargout{1} = 1;
           end
           return;
   end
 
   % check all the optional arguments
   if (mod(length(varargin),2) == 1)
       if (DISPLAY)
           disp ('ERROR: Invalid optional argument structure.');
       else
           varargout{1} = 2;
       end
       return;
   end

   % set the defaults
   samplingRate = 8000;
   sync = sin(pi/2*(0:31)');
   fftSize = 64;
   symbolsPerPacket = 14;
   cpLength = 8;
   
   % parse the optional arguments
   for i = 1:length(varargin)-1
       switch(varargin{i})
           case 'fs'
               samplingRate = varargin{i+1};
           case 'sync'
               sync = varargin{i+1};
           case 'fft_size'
               fftSize = varargin{i+1};
           case 'sym_per_packet'
               symbolsPerPacket = varargin{i+1};
           case 'cp_length'
               cpLength = varargin{i+1};
           otherwise
               if(DISPLAY)
                   disp('ERROR: Unrecognized optional argument.');
               else
                   varargout{1} = 2;
               end
               return;
       end
   end

   % calculate some other parameters we'll use
   % For our system we don't use all of the available spots.  
   % This formula gives us 24 carriers with a 64 point FFT.  
   % I'm not sure if it'll work for other FFT sizes.
   % Data will be in the format [0 D1 D2 ... D24 0 ... 0 D1* ... D24*]
   carriersPerSymbol = floor(0.4 * fftSize - 1);      
   bytesPerSymbol = carriersPerSymbol * bitsPerCarrier / 8;  
   bytesPerPacket = symbolsPerPacket * bytesPerSymbol;
   
   % read the input data (currently only supports ascii characters)
   fileID = fopen(in_file, 'r');
   if (fileID ~= -1)
      myData = fread(fileID, 'uchar');      
      fclose(fileID);
   else
      if (DISPLAY)
          disp ('ERROR: Could not open input file.');
      else
          varargout{1} = 3;
      end
      return;
   end

   % initialize the output with a quarter second of silence in beginning
   encodedData = zeros(samplingRate / 4,1);
   
   % Generate the time domain waveform in several steps
   %  (1) Figure out the total number of packets
   %  (2) For each packet...
   %     (2a) Add the sync code at the beginning
   %     (2b) Followed by channel estimation
   %     (2c) Followed by the symbol data
   %  (3) Optionally distort the time domain data with noise, freq.
   %      offset, channel effects, etc.
   %  (4) Write the data to a .wav file at the given sampling rate
   
   % (1) Number of packets needed to encode entire text file
   numEntirePackets = ceil(length(myData) / bytesPerPacket);  

   % If the length of myData is not an integer multiple of packet size, 
   % pad zeros to myData.
   if (length(myData) < numEntirePackets * bytesPerPacket)
      myData = [myData; ...
           zeros(numEntirePackets * bytesPerPacket - length(myData), 1)];
   end

   % (2) Encode each packet into the data stream
   for packetIndex = 1 : numEntirePackets
      % Put random gap inbetween packets
      minGap = 32;
      maxGap = 82;
      zeroGap = zeros(floor(rand() * (maxGap - minGap)) + minGap, 1);  
      encodedData = [encodedData; zeroGap];

      % (2a) Add the sync signal
      encodedData = [encodedData; sync];
    
      % (2b) Add the channel estimation
      % We transmit known data on each of the channels, where the known 
      % data is (1 + j*0).
      x = ones(carriersPerSymbol, 1);
      y = zeros(carriersPerSymbol, 1);

      % Equalization signal in frequency domain
      % Enforce complex-conjugate symmetry in order to guarantee real 
      % signals in the time domain.
      equalSignalInFreq = zeros(fftSize, 1); 
      equalSignalInFreq(2 : carriersPerSymbol+1) = x + j*y;                
      equalSignalInFreq(end-carriersPerSymbol+1 : end) = flipud(x - j*y);  
      
      % Convert to time domain and discard small imaginary part
      equalSignalInTime = real(ifft(equalSignalInFreq, fftSize));    
      % Add cyclic prefix (prepend the end of the sequence to the front)
      equalSignalInTime = ...
        [equalSignalInTime(end-cpLength+1:end); equalSignalInTime];  

      % Add the channel estimate to the encoded data
      encodedData = [encodedData; equalSignalInTime];

      % (2c) Generate OFDM data symbols based on modulation type
      % See Modulation HO for more details on how data is encoded
      dataIdx = 1+(packetIndex-1)*bytesPerPacket;
      for symbolIndex = 1 : symbolsPerPacket    
         % Real part of frequency domain carrier data
         x = zeros(carriersPerSymbol, 1);  
         % Imag part of frequency domain carrier data
         y = zeros(carriersPerSymbol, 1);  

         % get the bytes we need for this symbol
         symbolData = myData(dataIdx:dataIdx + bytesPerSymbol - 1);
         dataIdx = dataIdx + bytesPerSymbol;
         % convert it to a stream of bits
         symbolBits = [];
         for i = 1:length(symbolData)
            byte = dec2bin(symbolData(i));
            byte = byte - '0';
            symbolBits = [symbolBits; zeros(8 - length(byte),1); byte'];
         end
         
         % encode the bit stream into carriers
         carrierIdx = 1;
         for i = 1:bitsPerCarrier:length(symbolBits)
             bits = symbolBits(i:i+bitsPerCarrier-1);
             lookupIdx = 1;
             for k = 1:bitsPerCarrier
                 lookupIdx = lookupIdx + bits(k)*2^(bitsPerCarrier-k);
             end
             x(carrierIdx) = carrierLookup(lookupIdx, 1);
             y(carrierIdx) = carrierLookup(lookupIdx, 2);
             carrierIdx = carrierIdx + 1;
         end
         
         % Fill in zeros for unused space in the frequency domain and 
         % enforce complex-conjugate symmetry in order to guarantee real 
         % signals in the time domain.
         symbolInFreq = [0; x + j*y; ...
             zeros(fftSize-2*carriersPerSymbol-1, 1); flipud(x - j*y)];
         symbolInTime = real(ifft(symbolInFreq));
        
         % Add the cyclic prefix
         symbolInTime = [symbolInTime(end-cpLength+1:end); symbolInTime];

         % Append it to the encoded data
         encodedData = [encodedData; symbolInTime];
      end
   end

   % (3) Add optional channel effects, scaling, noise, etc.
   % -- TODO --
   % Make these effects options the user can set
   % -- TODO --
   % Scale for quantization to avoid overflow
   scaleFactor = 0.5; 
   encodedData = encodedData * scaleFactor;

   % Add AWGN noise to signal
   noisePower = 0.001;
   encodedData = encodedData + ...
       randn(length(encodedData),1) * noisePower * scaleFactor;

   % Add channel model
   h = [1 0.25 0.01];
   encodedData = conv(encodedData, h);

   % -- TODO --
   %freqOffsetPPM = 0;      % Frequency offset in parts per million 
   % Leave "freqOffsetPPM" to zero.  Actually, this offset feature
   % doesn't work properly.
   %freqOffset = freqOffsetPPM / 1e6; % Convert to fraction
   % resulting sampling frequency used to write wavefile
   %samplingRateWithFreqOffset = round(samplingRate * (1 + freqOffset)); 
   % -- TODO --
   
   % (4) Write to encoded data to a wavefile for playback to DSP board
   wavFileName = sprintf('%s.wav', out_file);
   wavwrite(encodedData, samplingRate, wavFileName);
   
   if(DISPLAY)
       disp('Data encoding complete.');
   else
       varargout{1} = 0;
   end
   return;
end