% Eric Bradley
% Image Compression Final Project
% Analog TV data to video/audio
% Note that this program has high memory requirements.  Recommend that you
% close all other programs and set the following switch to free up more
% virtual memory (works for XP):
% http://www.microsoft.com/technet/prodtechnol/exchange/guides/E2k3Perf_ScalGuide/e834e9c7-708c-43bf-b877-e14ae443ecbf.mspx?mfr=true
% The path structure is hard-coded.  
%    Bradley data is:  .\data\VCR_CH3_Video_3BLOCKS.dat
%    Beser data is:    .\data2\videoseg000.dat
%    Sub-directories  .\video  and  .\audio  must exist for creation of video and audio
 
clear;      % Free up all memory
close all;  % Close current figures
tic;        % Start the timer
diary on;   % Save console contents in case we crash (which we do sometimes with these large memory requirements
% Set current directory to where this file is
filename = mfilename('fullpath');
[pathstr, name, ext, ver] = fileparts(filename);
cd (pathstr);

% User defined booleans/variables
isEBData  = 1 ;  % 1=>Data from EWB NRL; 0=>Data from Dr. Beser
iter      = 1000 ;  % number of iterations to read in norig samples - processing will stop at EOF
doPlots   = 1 ;  % Boolean - display plots only on the first iteration?
saveTime  = 1 ;  % Boolean - don't do extraneous time consuming stuff?
doAudio   = 1 ;  % Boolean - should we extract the audio?
doVideo   = 1 ;  % Boolean - should we extract the luminance video?
doColor   = 1 ;  % Boolean - should we extract the color?
saveFiles = 1 ;  % Boolean - save extracted audio and video to disk?

% Data set specific inconsistencies
if (isEBData)
  fcaoffset = -0.0043e6;   % Audio peak is at 5.7457e6 in the data
  fifoffsetbeser = 0;
else
  fcaoffset = 0.0007e6;  % Audio peak is at 5.7507e6 in the data
  % In Beser's data:
  % The RF signal appears to be 0.06205 MHz less than expected.
  % The luminance carrier is at 6.18795 MHz in the RF instead of the expected 6.25
  % Color carrier 9.7676, exp 9.83
  % Audio carrier 10.687, exp 10.75
  % This is a good starting point, but the carrier freq appears to change
  % throughout the data - this code does not tolerant much change in
  % carrier (10s of Hz will cause severe degradation)
  fifoffsetbeser = -0.062016e6;  % THIS IS CRITICAL for proper demod
end 
% In Bradley's data:
% Audio carrier 4.4957, exp 5.75, diff 1.2543
% Color carrier 3.5795, exp 4.83, diff 1.2505
fifoffsetbradley = 1.25e6;

% Constants
lastDataSet = 0;  % Flag to indicate we're done reading
fcv = 1.25e6;  % Video Carrier Freq in IF
fcc = 3.579545e6 + fcv;  % Color Carrier Freq in IF
fca = 4.5e6 + fcv + fcaoffset;  % Audio Carrier Freq in IF
norig  = 1000000;  % Samples to read in
npl = 100       ;  % Num samples for a low density pretty plot
npm = 1000      ;  % Num samples for a medium density pretty plot
nph = 10000     ;  % Num samples for a high density pretty plot
fc = 60e6       ;  % RF Carrier Freq
Ba = 3*10e3     ;  % Bandwidth of audio (2 or 3 times 10e3 works well)
fcal = fca - Ba ;
fcah = fca + Ba ;
fccl = 2.1e6 + fcv ;  % BPF cutoff
fcch = 4.1e6 + fcv ;  % BPF cutoff
Bc = 6.0e6      ;  % Bandwidth of entire signal
Bv = 4.2e6      ;  % Bandwidth of video
Bci = 1.5e6  ;  % Bandwidth of I color  
Bcq = 0.5e6  ;  % Bandwidth of Q color  
butterOrder = 5 ;  % if we go too big (>5) we start to lose the audio 

% In order to get Bradley's spectrum to look like Beser's (to ease
% computation) we must upconvert it
fup = 5e6 + fifoffsetbradley + fifoffsetbeser;    % upconversion frequency to get to IF 
% IF frequency band
fifl = 5e6 + fifoffsetbeser;
fifh = fifl + Bc;
% Analog sample frequency is normalized for each data set to be the same as
% Beser data for simplicity in data manipulation
fs = 32e6;  % Analog sampling freq
Ts = 1/fs     ;  % Analog sampling period
% Audio sample rates
fsa = 44100   ;  % Sampling frequency of audio signal
Tsa = 1/fsa   ;
% Field parameters
fld = 0;  % cummulative count of number of fields
% left over video signal (partial field) from previous iteration
leftOverV = [];  
leftOverCI = [];  
leftOverCQ = [];  
% audio signal initialization
audio = [];

if (isEBData)
  % Read the raw data in - IF 0 - 6 MHz
  % Binary file format:  3 doubles - dt, offset, gain
  %                      A[0] = I[0] - signed 16-bit integer
  %                      A[1] = Q[0] - signed 16-bit integer
  %                      A[2] = I[1] - signed 16-bit integer
  %                      A[3] = Q[1] - signed 16-bit integer
  %                      ...          
  inputFileName = 'VCR_CH3_Video_3BLOCKS.dat'; % Each block is a second in length!  :(
  %inputFileName = 'VCR_CH3_NO_Video_01.dat';
  %inputFileName = 'VCR_CH3_Video_01.dat';
  fid = fopen(['data\' inputFileName], 'r', 'ieee-be');  % Big Endian
  params = fread(fid,  3, 'double', 0); 
  Tsbradley = params(1);
  fsbradley = 1/Tsbradley;
  gain = params(3);
else
  inputFileName = 'videoseg000.dat';
  fid = fopen(['data2\' inputFileName], 'r', 'ieee-be'); % Big Endian
  % read the 512 byte header, and never look at it again
  [header, cnthead] = fread(fid, 256, 'int16=>int16');
end
% Baseband fs and Ts
fsb = 2*Bc;
Tsb = 1/fsb;

% nth order Butterworth filters 
%   needn't do this every iteration, so outside loop
% IF to baseband band pass for Beser data
[b_if_bpf, a_if_bpf] = butter(butterOrder, [fifl fifh]/(fs/2));
% IF to baseband low pass for Beser data
[b_if_lpf, a_if_lpf] = butter(butterOrder, Bc/(fs/2), 'low');
% BPF Audio 
[b_aud_bpf, a_aud_bpf] = butter(butterOrder, [fcal fcah]/(fsb/2));   
% LPF Video Luminance - 4.2 MHz above luminance center (1.25 MHz)
[b_vid_lpf, a_vid_lpf] = butter(butterOrder, (Bv)/(fsb/2), 'low');
% BPF Color - 2.1-4.1 MHz above luminance center (1.25 MHz)
[b_col_bpf, a_col_bpf] = butter(butterOrder, [fccl fcch]/(fsb/2));
% LPF Color I - 1.5 MHz
[b_i_lpf, a_i_lpf] = butter(butterOrder, Bci/(fsb/2), 'low');
% LPF Color Q - 0.5 MHz
[b_q_lpf, a_q_lpf] = butter(butterOrder, Bcq/(fsb/2), 'low');

for loop=1:iter
% while(cntdata>0)
%     [data,cntdata]=fread(fid1,1047552*4,'int16=>int16');
if (isEBData)
  [yrf2 nread] = fread(fid, norig*2, 'int16');
  if nread ~= norig*2
    lastDataSet=1; 
    if nread == 0, break; end;
  end
  yrf2 = yrf2';
  % Data was sampled in IQ (QAM) modulation.  We want both "channels", so combine
  I = yrf2(1:2:nread);
  Q = yrf2(2:2:nread);
  yrf2 = sqrt(I.^2 + Q.^2);
  % Note that the data has a DC bias.  We will not compensate for that here
  % as the frequency shifting will remove the bias
  
  % upsample so both data sets have the same sample period
  yrf = resample(yrf2, fs, fsbradley);  
  nread = norig * fs/fsbradley;  % nread changes due to resample
else
  [yrf nread] = fread(fid, norig, 'int16');
  if nread ~= norig
    lastDataSet=1; 
    if nread == 0, break; end;
  end
  yrf = yrf';
end

nb = nread * (fsb/fs);  % # samples @ baseband

% Plot RF data
t = Ts * (0:nread-1);
f = fs/2 * linspace(0,1,nread/2-1);
Yrf = fft(yrf);
if (doPlots) 
  figure('Name', 'RF Signal');
  subplot(3,1,1), plot(t(1:npl), yrf(1:npl), 'x-'), title('RF Waveform');
  subplot(3,1,2), plot(t(1:nph), yrf(1:nph)), title('RF Waveform');
  subplot(3,1,3), plot(f, log10(abs(Yrf(2:nread/2)))), title('RF Freq Response');
end;  

% BPF at IF, Downconvert to baseband, and LPF
if (isEBData)  
  % Shift it up to 5-11 MHz IF where the other data set resides, 
  % and from that point the processing is the same.
  yrfup = yrf .* cos(2*pi*fup*t);   % Upconvert to IF of 5-11 MHz
else
  yrfup = yrf;  % Beser data doesn't need upconversion
end
% Band pass filter to eliminate other harmonics when mixing
yrffilt = filter(b_if_bpf, a_if_bpf, yrfup);
if not(saveTime)
  strLabel = sprintf('IF to Baseband %dth Order Butterworth Filter', butterOrder);
  figure('Name', strLabel), freqz(b_if_bpf, a_if_bpf, norig, fs), title('Band pass before mixing'); 
end
% Data has been upconverted into the 5-11 MHz IF band
% Mix it to baseband
yif = yrffilt .* cos(2*pi*fifl*t);   % Mix with carrier freq  
% Low pass filter to remove unwanted harmonics
yfilt = filter(b_if_lpf, a_if_lpf, yif);
if not(saveTime)
  strLabel = sprintf('IF to Baseband %dth Order Butterworth Filter', butterOrder);
  figure('Name', strLabel), freqz(b_if_lpf, a_if_lpf, nread, fs), title('Low pass after mixing');
end;
% FFT all the signals  
Yrffilt = fft(yrffilt);
Yif = fft(yif);
Yfilt = fft(yfilt);
f = fs/2 * linspace(0,1,nread/2-1);
if (doPlots) 
  figure('Name', 'IF Signal');
  subplot(3,1,1), plot(f, log10(abs(Yrffilt(2:nread/2)))), title('Filtered IF Freq Response');
  subplot(3,1,2), plot(f, log10(abs(Yif(2:nread/2)))), title('Mixed IF Freq Response');
  subplot(3,1,3), plot(f, log10(abs(Yfilt(2:nread/2)))), title('LPF Baseband Freq Response');
end

% Resample the data at the Nyquist rate at baseband to make future calculations more efficient
y = resample(yfilt, fsb, fs);
Y = fft(y);
f = fsb/2 * linspace(0,1,nb/2-1);
t = Tsb * (0:nb-1);
if (doPlots) 
  figure('Name', 'Baseband Signal');
  subplot(2,1,1), plot(f, log10(abs(Y(2:nb/2)))), title('Baseband Freq Response');
  subplot(2,1,2), plot(t(1:npm), y(1:npm)), title('Baseband Waveform');
end

% At baseband, v(t) is B&W video, a(t) is audio, c(t) is color video

% Extract Audio 
if (doAudio)
  % BPF Audio 
  ya = filter(b_aud_bpf, a_aud_bpf, y);
  Ya = fft(ya);

  if (doPlots) 
    if not(saveTime)
      strLabel = sprintf('Audio Signal %dth Order Butterworth Filter', butterOrder);
      figure('Name', strLabel);
      freqz(b_aud_bpf, a_aud_bpf, nb, fsb);
    end;
    figure('Name', 'Modulated Audio Signal');
    subplot(3,1,1), plot(t(100:100+npl), ya(100:100+npl)), title('Modulated Audio Waveform');
    subplot(3,1,2), plot(t(100:100+npm), ya(100:100+npm)), title('Modulated Audio Waveform');
    subplot(3,1,3), plot(f, log10(abs(Ya(2:nb/2)))), title('Modulated Audio Freq Response');
  end;  

  % Demodulate audio:  FM (frequency modulation)
  aa = demod(ya, fca, fsb, 'fm');
  a = resample(aa, fsa, fsb);  % downconvert to audio-type frequencies
  a(1:10) = 0;   % there is always a blip at the beginning.  zero it out
  na = numel(a);
  ta = Tsa * (0:na-1);
  fa = fsa/2 * linspace(0,1,na/2-1);
  A = fft(a);
  if (doPlots)
    figure('Name', 'Audio Signal');
    subplot(2,1,1), plot(ta, a), title('Audio Waveform');
    subplot(2,1,2), plot(fa, log10(abs(A(2:na/2)))), title('Audio Freq Response');
  end;
  audio = [audio a];   % store audio in array

end % if doAudio

% Extract Video
if (doVideo)
  % Demodulate video:  VSB + C (vesidual sibe band + carrier)
  % We have already downconverted the signal to remove the carrier frequency
  yv = demod(y, fcv, fsb, 'amssb');
  % Beser data is inverted...sometimes
  if not(isEBData) && min(yv) < -50 
    yv = -yv;  
  end
  Yv = fft(yv);

  % LPF to get v(t)
  % Note: A comb filter would be better for v(t), but do LPF for now
  v = filter(b_vid_lpf, a_vid_lpf, yv);
  V = fft(v);
  if (doPlots)
    if not(saveTime)
      strLabel = sprintf('Video Signal %dth Order Butterworth Filter', butterOrder);
      figure('Name', strLabel);
      freqz(b_vid_lpf, a_vid_lpf, nb, fsb);
    end;
    figure('Name', 'Luminance Video Signal');
    subplot(3,1,1), plot(t(1:nph), v(1:nph)), title('Luminance Video Waveform');
    subplot(3,1,2), plot(t, v), title('Luminance Video Waveform');  
    subplot(3,1,3), plot(f, log10(abs(V(2:nb/2)))), title('Luminance Video Freq Response');
    figure('Name', 'Luminance Video Waveform'), strips(v(1:100000),5000*Tsb,fsb), title('Luminance Video Waveform');
  end

  % Filter and demod color
  if (doColor)
    c = filter(b_col_bpf, a_col_bpf, y); 
    C = fft(c);

    if (doPlots) 
      if not(saveTime)
        strLabel = sprintf('Color Signal %dth Order Butterworth Filter', butterOrder);
        figure('Name', strLabel);
        freqz(b_col_bpf, a_col_bpf, nb, fsb);
      end;
      figure('Name', 'Modulated Color Signal');
      subplot(3,1,1), plot(t(1:npm), c(1:npm)), title('Modulated Color Waveform');
      subplot(3,1,2), plot(t(1:nph), c(1:nph)), title('Modulated Color Waveform');
      subplot(3,1,3), plot(f, log10(abs(C(2:nb/2)))), title('Modulated Color Freq Response');
    end;  

    % Color demod - QAM
    [ciunfilt, cqunfilt] = demod(c, fcc, fsb,'qam');
    ci = filter(b_i_lpf, a_i_lpf, ciunfilt);
    cq = filter(b_q_lpf, a_q_lpf, cqunfilt);  
    CI = fft(ci);
    CQ = fft(cq);
    if (doPlots)
      figure('Name', 'Color Video Signal');
      subplot(4,1,1), plot(t, ci), title('Color I Video Waveform');  
      subplot(4,1,2), plot(t, cq), title('Color Q Video Waveform');  
      subplot(4,1,3), plot(f, log10(abs(CI(2:nb/2)))), title('Color I Video Freq Response');
      subplot(4,1,4), plot(f, log10(abs(CQ(2:nb/2)))), title('Color Q Video Freq Response');
      figure('Name', 'Color I Video Waveform'), strips(ci(1:100000),5000*Tsb,fsb), title('Color I Video Waveform');
      figure('Name', 'Color Q Video Waveform'), strips(cq(1:100000),5000*Tsb,fsb), title('Color Q Video Waveform');
    end
  end % if doColor

  % 30 frames/second (33.3 ms/frame) (2 interlaced fields/frame), 
  % 60 fields/second (16.7 ms/field) (actually 59.95, but we'll ignore that for now)
  % 525 lines/frame, 262.5 lines/field (1st field ends at 1/2, 2nd field starts at 1/2
  % 15750 lines/second, 6.349e-5 (63.49 us) sec/line
  % 63.5 us line time: 10 us horizontal retrace, 53.5 us active line time
  % @ 4:3 ratio, 700 pels/line in 53.5 us = 7.643e-8 (76.43 ns) sec/pel 
  % Vertical sync pulse: 27 us
  
  % Video stats
  if isEBData
    nlines = 262;   % lines per field WARNING for now throw away the extra 1/2 line
  else % since we have not vsync for Beser data, keep the vsync as if its part of the field
    nlines = 262+21; % 262.5 + 20.5 (tvblank/(tactline+thretrace))
  end
  nlinesframe = 2*nlines;
  npels = 700;    % pels per line
  nvsyncpul = 6;  % number of vertical resync pulses
  % Video signal timing
  tactline = 53.5e-6;
  thretrace = 10e-6;
  tpel = tactline / npels;
  tvpulse = 27e-6; 
  tvpulseh = tvpulse * 1.5;  % allow some slop
  tvpulsel = tvpulse * 0.5;
  tfield = 1 / 59.94;
  tvblank = 1.3e-3;  % approximate time of entire vertical blanking period (tfield*0.07)
  % Video Levels
  lhsync = 100;
  lhsyncstartmin = 75;%lsync - 10;
  lhsyncend = 68;  % if this is too close to 75, you will pick up the front of the back porch instead of the end
  lvsync = 85;
  lmax = lhsync;
  lblack = 75;
  lwhite = 12.5;
  % Color conversion
  yiq2rgb = [1.0  0.956  0.620;
             1.0 -0.272 -0.647;
             1.0 -1.108  1.700 ];
  % Normalize to 0-100 - assume it is all positive to start
  % find the max of the non-end points - these can still be outliers for some reason
  v = v * lmax/max(v(100:numel(v)-100));  % this could give inconsistencies between iterations, but seems to be OK
  % Trim down those that are out of range - the first and last 5 often are
  out = find(v > 100);
  v(out) = 100;
  out = find(v < 0);
  v(out) = 0;

  % Add the left over from last time
  v = [leftOverV v];
  if doColor
    ci = [leftOverCI ci];
    cq = [leftOverCQ cq];
  end
  tv = Tsb * (0:numel(v)-1);
  if (doPlots)
    figure('Name', 'Normalized Luminance Video Signal');
    subplot(2, 1, 1), plot(tv, v), title('Normalized Luminance Video Waveform');  
    subplot(2, 1, 2), plot(tv(1:nph), v(1:nph)), title('Normalized Luminance Video Waveform');  
  end

  
  i = 1;    % i is current index into video signal   
  while 1  % loop is terminated when a full field is not found
    % Rewind a full vertical blanking duration to find next vertical sync
    % this is extra conservative since we shouldn't be at the end of the
    % vert blanking
    i = i - tvblank / Tsb;  
    if i <= 0, i = 1; end;
    % Find a vertical sync pulse sequence
    if isEBData  % WARNING no vert sync for Beser data
      foundVert = 0;
      while foundVert == 0
        for j=1:nvsyncpul   % count the number of pulses
          i2 = i - 1 + find(v(i:numel(v)) > lvsync, 1);  % find next rising edge over vsync threshold
          i = i2 - 1 + find(v(i2:numel(v)) < lvsync, 1); % find next falling edge
          tp = (i - i2) * Tsb;  % measure the pulse width
          if (tp > tvpulseh || tp < tvpulsel)    % not the pulse
            break;   % start counting over
          end
        end % for
        if j == nvsyncpul, foundVert = 1; end;
      end
    end
    
    % Check to see if there is enough time left to get a full field (plus
    %   some fudge factor
    if (numel(v)-i < 1.1 * tfield / Tsb)
      % Save the left over video signal (partial field) plus some extra on
      %   the front end for next time
      leftOverV = v(i-tvblank/Tsb:numel(v));
      if (doColor)
        leftOverCI = ci(i-tvblank/Tsb:numel(ci));
        leftOverCQ = cq(i-tvblank/Tsb:numel(cq));
      end
      break;
    end;
    fld = fld + 1;

    for l=1:nlines
      % Find a horizontal retrace sync pulse
      iorig = i;
      i = i - 1 + find (v(i:numel(v)) > lhsyncstartmin, 1);
      
      % wait for the entire delay duration or find the edge
      % we do both because the delay is too much sometimes
      % this isn't perfect.  we still end up grabbing a little too
      % much and get into the horizontal sync, but maybe that's how it works
      i1 = i + thretrace/Tsb;
      i2 = i - 1 + find(v(i:numel(v)) < lhsyncend, 1);
      i = min(i1, i2);
      if isempty(i)
        i = iorig + (nlines-l) * tactline/Tsb;  % skip the rest of the field
        disp(['Warning: Line ' int2str(l) ' of Field ' int2str(fld) ' was unreadable.  Skipping field.']);
        break;  % if the data is really noisy, give up on the field
      end      
      if i + tactline/Tsb > numel(v)
        i = iorig + (nlines-l) * tactline/Tsb;  % skip the rest of the field
        disp(['Warning: Line ' int2str(l) ' of Field ' int2str(fld) ' was too long.  Skipping field.']);
        break;  % if the data is really noisy, give up on the field
      end      

      % Extract the duration of the line
      lineraw = v(i : i+(tactline/Tsb));
      % Resample to the right number of pels with nearest neighbor interpolation
      line = resample(lineraw, npels, numel(lineraw), 0);
      % Invert and Normalize the color - 0=black, 255=white
      line = (lblack-line)*255/(lblack-lwhite);
      % Build up the field image
      field(l, :, fld) = uint8(line);  % this is crucial to be uint8 to save memory - large contiguous storage location
      
      if (doColor)
        lineiraw = ci(i : i+(tactline/Tsb));
        lineqraw = cq(i : i+(tactline/Tsb));
        linei = resample(lineiraw, npels, numel(lineiraw), 0);
        lineq = resample(lineqraw, npels, numel(lineqraw), 0);
        fieldi(l, :) = linei;  % these need to be doubles - need neg #s
        fieldq(l, :) = lineq;
      end

      % Increment for the next line read attempt
      i = i + tactline / Tsb;
    end % l

    if doPlots && fld == 2
      frame(1:2:nlinesframe,:) = field(:,:,1);
      frame(2:2:nlinesframe,:) = field(:,:,2);
      figure('Name', 'Luminance Fields and Frame');
      subplot(3,1,1), imshow(field(:,:,1)), title('Luminance Field 1');
      subplot(3,1,2), imshow(field(:,:,2)), title('Luminance Field 2');
      subplot(3,1,3), imshow(frame), title('Luminance Frame');
    end

    % Translate field to RGB
    if (doColor)
      for yi=1:nlines
        for xi=1:npels
          fieldRGB(yi,xi,:,fld)=uint8(yiq2rgb*[double(field(yi,xi,fld));fieldi(yi,xi);fieldq(yi,xi)]);
        end
      end
      if doPlots && fld == 2
        % I & Q are bipolar, so scale so they are viewable
        figure('Name', 'I, Q, and RGB Fields');
        subplot(3,1,1), imshow(uint8((fieldi+256)/2)), title('Normalized I Field');
        subplot(3,1,2), imshow(uint8((fieldq+256)/2)), title('Normalized Q Field');
        subplot(3,1,3), imshow(fieldRGB(:,:,:,1)), title('RGB Field');
      end
    end

  end % while 1

end % if doVideo

% print out loop iteration so we know where we are
disp(['Iteration ' int2str(loop) ' of ' int2str(iter) ' finished.  (Field # ' int2str(fld) ')'])  
toc;
doPlots = 0;   % Only do it for the first iteration

if lastDataSet == 1, break; end;

end  % end for loop

fclose(fid);


if (doAudio)
  % Combine all audio segments
  audio=audio/max(abs(audio));    % normalize it
  if saveFiles
    audFileName = ['audio/' inputFileName(1:length(inputFileName)-4) '_' datestr(now,'dd-mmm-yy-HH-MM-SS') '.wav'];
    wavwrite(audio, fsa, audFileName)
  end % if saveFiles
  sound(audio, fsa)
end % if doAudio

if (doVideo)
  cmapBW = colormap(gray(256));
  for i=1:floor(fld/2)  % this will throw away an odd numbered last field
    % Deinterlace - each field contains some negative bits, but those are only off screen as expected
    frame(1:2:nlinesframe,:) = field(:,:,i*2-1);
    frame(2:2:nlinesframe,:) = field(:,:,i*2);
 
    if saveFiles || i <= 50
      % Create frame format required for movie2avi
      mov(i) = im2frame(frame(:,:), cmapBW);
    end
    
    if i <= 50  % MATLAB movie crashes if too big
      % MATLAB movie format is flipped vertically
      movfl(i) = mov(i);
      movfl(i).cdata = flipud(mov(i).cdata);
    end
    
    if doColor
      frameRGB(1:2:nlinesframe,:,:) = fieldRGB(:,:,:,i*2-1);
      frameRGB(2:2:nlinesframe,:,:) = fieldRGB(:,:,:,i*2);
      movRGB(i) = im2frame(frameRGB(:,:,:));
    end
  end
  
  % save movie
  if saveFiles
    vidFileName = ['video/' inputFileName(1:length(inputFileName)-4) '_' datestr(now,'dd-mmm-yy-HH-MM-SS') '.avi'];
    % Save uncompressed video - 
    %  if this becomes a memory problem, use colormap(gray(64)) on im2frame which is 64
    %  instead of 256 entries, and leave off the compression param on
    %  movie2avi to get decent compression
    movie2avi(mov, vidFileName, 'fps', 30, 'compression', 'None');

    if doColor
      vidRGBFileName = ['video/' inputFileName(1:length(inputFileName)-4) '_RGB_' datestr(now,'dd-mmm-yy-HH-MM-SS') '.avi'];
      % Save uncompressed video - 
      %  if this becomes a memory problem, use colormap(jet(64)) on im2frame which is 64
      %  instead of 256 entries, and leave off the compression param on
      %  movie2avi to get decent compression
      movie2avi(movRGB, vidRGBFileName, 'fps', 30, 'compression', 'None');
    end
  end % if saveFiles  
  
  % Show movie - this fails sometimes, so put it at the end
  fig = figure('position', [200 200 npels nlinesframe], 'Name', 'Video Sample - first 50 Frames');
  % For some reason movie() craps out with more than 50 frames - or sometimes less
  % See bug: http://www.mathworks.com/support/bugreports/details.html?rp=319426
  movie(fig, movfl(1:min(50,floor(fld/2))), -5, 30);  % play movie forwards & back 5 times at 30 FPS

end  % if doVideo

toc;
return

