% 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