From fa03578f04d51367a4fea1306037ae883ecc6213 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Mon, 12 Mar 2018 20:58:07 +0000 Subject: [PATCH] restored ofdm Octave/C code from last working rev 3387, tofdm passing ok git-svn-id: https://svn.code.sf.net/p/freetel/code@3408 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/ofdm_lib.m | 42 +- codec2-dev/octave/tofdm.m | 534 +++++------------------ codec2-dev/src/codec2_ofdm.h | 4 +- codec2-dev/src/ofdm.c | 765 +++++---------------------------- codec2-dev/src/ofdm_internal.h | 105 +++-- codec2-dev/unittest/tofdm.c | 225 ++++------ 6 files changed, 376 insertions(+), 1299 deletions(-) diff --git a/codec2-dev/octave/ofdm_lib.m b/codec2-dev/octave/ofdm_lib.m index b6839993..32dd02f1 100644 --- a/codec2-dev/octave/ofdm_lib.m +++ b/codec2-dev/octave/ofdm_lib.m @@ -63,7 +63,7 @@ endfunction + more suitable for real time implementation #} -function [t_est foff_est ratio] = coarse_sync(states, rx, rate_fs_pilot_samples) +function [t_est foff_est] = coarse_sync(states, rx, rate_fs_pilot_samples) Nsamperframe = states.Nsamperframe; Fs = states.Fs; Npsam = length(rate_fs_pilot_samples); verbose = states.verbose; @@ -71,18 +71,13 @@ function [t_est foff_est ratio] = coarse_sync(states, rx, rate_fs_pilot_samples) Ncorr = length(rx) - (Nsamperframe+Npsam) + 1; assert(Ncorr > 0); corr = zeros(1,Ncorr); - mag_max = 0; for i=1:Ncorr corr(i) = abs(rx(i:i+Npsam-1) * rate_fs_pilot_samples'); corr(i) += abs(rx(i+Nsamperframe:i+Nsamperframe+Npsam-1) * rate_fs_pilot_samples'); - mag_l = max( abs(rx(i)) , abs(rx(i+Nsamperframe)) ); - mag_max = max(mag_max, mag_l); end [mx t_est] = max(abs(corr)); - ratio = (mx + 1e-9) / (mag_max + 1e-9); - C = abs(fft(rx(t_est:t_est+Npsam-1) .* conj(rate_fs_pilot_samples), Fs)); C += abs(fft(rx(t_est+Nsamperframe:t_est+Nsamperframe+Npsam-1) .* conj(rate_fs_pilot_samples), Fs)); @@ -96,8 +91,8 @@ function [t_est foff_est ratio] = coarse_sync(states, rx, rate_fs_pilot_samples) foff_est = foff_est_neg - fmax - 1; end - %printf("t_est: %d\n", t_est); if verbose > 1 + %printf("t_est: %d\n", t_est); figure(7); clf; plot(abs(corr)) figure(8) @@ -173,7 +168,6 @@ function states = ofdm_init(bps, Rs, Tcp, Ns, Nc) states.sample_point = states.timing_est = 1; states.nin = states.Nsamperframe; - states.sync = 0; % generate OFDM pilot symbol, used for timing and freq offset est rate_fs_pilot_samples = states.pilots * W/states.M; @@ -434,38 +428,6 @@ function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod(states, states.foff_est_hz = foff_est_hz; endfunction -function [rx_bits states aphase_est_pilot_log rx_np rx_amp] = ofdm_demod2(states, rxbuf_in) - ofdm_load_const; - [rx_bits states aphase_est_pilot_log rx_np rx_amp] =ofdm_demod(states,rxbuf_in); - - - [t_est foff_est ratio] = coarse_sync(states,states.rxbuf,rate_fs_pilot_samples); - t_est = mod(t_est,Nsamperframe) - (M+Ncp); - printf("t_est is %d, foff_est is %d, ofdm_foff_est is %f, ratio %f\n",t_est,foff_est,states.foff_est_hz,ratio) - %If out of sync, try a coarse sync on RX buffer - if states.sync <= 0 - %[t_est foff_est ratio] = coarse_sync(states,states.rxbuf,rate_fs_pilot_samples); - %Correct to be within one frame and at head of frame instead of after first symbol - if(ratio > .5) - if(t_est > 12) - states.nin = t_est; - end - if(abs(foff_est - states.foff_est_hz) > 3) - states.foff_est_hz = foff_est; - end - printf("Syncing\n") - states.sync = 3; - end - else - if(ratio < .5) - states.sync = states.sync-1; - end - if(abs(foff_est - states.foff_est_hz) > 3) - states.foff_est_hz = foff_est; - end - end -endfunction - % Save test bits frame to a text file in the form of a C array % diff --git a/codec2-dev/octave/tofdm.m b/codec2-dev/octave/tofdm.m index 0d64e22a..73b99fea 100644 --- a/codec2-dev/octave/tofdm.m +++ b/codec2-dev/octave/tofdm.m @@ -5,461 +5,155 @@ % ------------------------------------------------------------------ -load rand_state_vec.mat -%rand("state",rand_state) - -function rx = ebno_awgn_channel(tx, Fs, Rb, EbNodB) - %rand("seed",2) - %randn("seed",3); - EbNo = 10^(EbNodB/10); - variance = Fs / (Rb*EbNo); - noise = sqrt(variance)*randn(1,length(tx)); - avg = sum(abs(tx))/length(tx); - rx = noise.*avg + tx; - %rx = rx * .10; -end +Nframes = 30; +sample_clock_offset_ppm = 100; +foff_hz = 0.5; -function rx = snr_awgn_channel(tx, Fs, BwS, snrdb) - snr = 10^(snrdb/10); - %signal power - psig = sum(abs(tx))^2 / length(tx); - %noise power in interesting band - pnoiseb = (1/snr) * psig; - %noise power in 1hz band - pnoise0 = pnoiseb/BwS; - %noise power over sample rate bandwidth - pnoiset = pnoise0 * (Fs/2); - - noise = sqrt(pnoiset)*randn(1,length(tx)); - - rx = noise + tx; - %rx = rx * .10; -end +more off; format; +ofdm_lib; +autotest; -function nums = im_re_interleave(nim) - nums = zeros(1,length(nim)*2); - nums(1:2:length(nums)) = real(nim); - nums(2:2:length(nums)) = imag(nim); +% --------------------------------------------------------------------- +% Run Octave version +% --------------------------------------------------------------------- -end +Ts = 0.018; Tcp = 0.002; Rs = 1/Ts; bps = 2; Nc = 16; Ns = 8; +states = ofdm_init(bps, Rs, Tcp, Ns, Nc); +ofdm_load_const; -% Trim two input arrays down to the shortest common length -function [a b] = trim1(ain, bin) - com = min(length(ain),length(bin)); - a = ain(1:com); - b = bin(1:com); -endfunction - -function [a b] = trim2(ain, bin) - com = min( size(ain)(1), size(bin)(1) ); - a = ain(1:com,:); - b = bin(1:com,:); -endfunction - -function pass = run_ofdm_test(Nframes,sample_clock_offset_ppm,foff_hz,EbNodB = 100,nheadsamp = 0) - rand("seed",2) - randn("seed",2) - %Nframes = 30; - %sample_clock_offset_ppm = 100; - %foff_hz = 5; - - more off; format; - ofdm_lib; - autotest; - - % --------------------------------------------------------------------- - % Run Octave version - % --------------------------------------------------------------------- - - Ts = 0.018; Tcp = 0.002; Rs = 1/Ts; bps = 2; Nc = 16; Ns = 8; - Rb_real = 700; - states = ofdm_init(bps, Rs, Tcp, Ns, Nc); - ofdm_load_const; - - tx_bits = round(rand(1,Nbitsperframe)); - - % Run tx loop - - tx_bits_log = []; tx_log = []; - for f=1:Nframes - tx_bits_log = [tx_bits_log tx_bits]; - tx_log = [tx_log ofdm_mod(states, tx_bits)]; - end +rand('seed',1); +tx_bits = round(rand(1,Nbitsperframe)); - % Channel simulation ---------------------------------------------- +% Run tx loop - rx_log = sample_clock_offset(tx_log, sample_clock_offset_ppm); - rx_log = freq_shift(rx_log, foff_hz, Fs); - #rx_log(1:nheadsamp) = zeros(nheadsamp,1); - rx_log_l = length(rx_log) - rx_log = [zeros(1,nheadsamp), rx_log]; - rx_log = rx_log(1:rx_log_l); - rx_log = ebno_awgn_channel(rx_log, Fs, Rb_real, EbNodB); +tx_bits_log = []; tx_log = []; +for f=1:Nframes + tx_bits_log = [tx_bits_log tx_bits]; + tx_log = [tx_log ofdm_mod(states, tx_bits)]; +end - rx_vec = fopen("tofdm_rx_vec","wb+"); - fwrite(rx_vec,im_re_interleave(rx_log),"single"); - fclose(rx_vec); +% Channel simulation ---------------------------------------------- - % Rx --------------------------------------------------------------- +rx_log = sample_clock_offset(tx_log, sample_clock_offset_ppm); +rx_log = freq_shift(rx_log, foff_hz, Fs); - % Init rx with ideal timing so we can test with timing estimation disabled +% Rx --------------------------------------------------------------- - Nsam = length(rx_log); - prx = 1; - nin = Nsamperframe+2*(M+Ncp); - %states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx_log(prx:nin); - %prx += nin; +% Init rx with ideal timing so we can test with timing estimation disabled - rxbuf_log = []; rxbuf_in_log = []; rx_sym_log = []; foff_hz_log = []; - timing_est_log = []; sample_point_log = []; - phase_est_pilot_log = []; rx_amp_log = []; - rx_np_log = []; rx_bits_log = []; +Nsam = length(rx_log); +prx = 1; +nin = Nsamperframe+2*(M+Ncp); +states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx_log(prx:nin); +prx += nin; - states.timing_en = 1; - states.foff_est_en = 1; - states.phase_est_en = 1; +rxbuf_log = []; rxbuf_in_log = []; rx_sym_log = []; foff_hz_log = []; +timing_est_log = []; sample_point_log = []; +phase_est_pilot_log = []; rx_amp_log = []; +rx_np_log = []; rx_bits_log = []; - if states.timing_en == 0 - % manually set ideal timing instant - states.sample_point = Ncp; - end +states.timing_en = 1; +states.foff_est_en = 1; +states.phase_est_en = 1; - have_frame = true; - while have_frame - % insert samples at end of buffer, set to zero if no samples - % available to disable phase estimation on future pilots on last - % frame of simulation - - nin = states.nin; - lnew = min(Nsam-prx+1,nin); - rxbuf_in = zeros(1,nin); - %printf("nin: %d prx: %d lnew: %d\n", nin, prx, lnew); - if lnew - rxbuf_in(1:lnew) = rx_log(prx:prx+lnew-1); - end - prx += nin; - if prx > Nsam - have_frame = false; - end - - [rx_bits states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod2(states, rxbuf_in); - % log some states for comparison to C - - rxbuf_in_log = [rxbuf_in_log rxbuf_in]; - rxbuf_log = [rxbuf_log states.rxbuf]; - rx_sym_log = [rx_sym_log; states.rx_sym]; - phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log]; - rx_amp_log = [rx_amp_log arx_amp]; - foff_hz_log = [foff_hz_log; states.foff_est_hz]; - timing_est_log = [timing_est_log; states.timing_est]; - sample_point_log = [sample_point_log; states.sample_point]; - rx_np_log = [rx_np_log arx_np]; - rx_bits_log = [rx_bits_log rx_bits]; - - end - - %states.foff_est_hz = 10; - % for f=1:Nframes +if states.timing_en == 0 + % manually set ideal timing instant + states.sample_point = Ncp; +end - % % insert samples at end of buffer, set to zero if no samples - % % available to disable phase estimation on future pilots on last - % % frame of simulation - - % nin = states.nin; - % lnew = min(Nsam-prx+1,nin); - % rxbuf_in = zeros(1,nin); - % %printf("nin: %d prx: %d lnew: %d\n", nin, prx, lnew); - % if lnew - % rxbuf_in(1:lnew) = rx_log(prx:prx+lnew-1); - % end - % prx += lnew; - - % [rx_bits states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod2(states, rxbuf_in); - - % % log some states for comparison to C - - % rxbuf_in_log = [rxbuf_in_log rxbuf_in]; - % rxbuf_log = [rxbuf_log states.rxbuf]; - % rx_sym_log = [rx_sym_log; states.rx_sym]; - % phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log]; - % rx_amp_log = [rx_amp_log arx_amp]; - % foff_hz_log = [foff_hz_log; states.foff_est_hz]; - % timing_est_log = [timing_est_log; states.timing_est]; - % sample_point_log = [sample_point_log; states.sample_point]; - % rx_np_log = [rx_np_log arx_np]; - % rx_bits_log = [rx_bits_log rx_bits]; - - % end - - - % --------------------------------------------------------------------- - % Run C version and plot Octave and C states and differences - % --------------------------------------------------------------------- - - % Override default path by setting path_to_tofdm = "/your/path/to/tofdm" - - if exist("path_to_tofdm", "var") == 0 - path_to_tofdm = "../build_linux/unittest/tofdm"; +for f=1:Nframes + + % insert samples at end of buffer, set to zero if no samples + % available to disable phase estimation on future pilots on last + % frame of simulation + + nin = states.nin; + lnew = min(Nsam-prx+1,nin); + rxbuf_in = zeros(1,nin); + %printf("nin: %d prx: %d lnew: %d\n", nin, prx, lnew); + if lnew + rxbuf_in(1:lnew) = rx_log(prx:prx+lnew-1); end - system(path_to_tofdm); - - load tofdm_out.txt; - % Generated with modem probe thing - load ofdm_test.txt; - - fg = 1; - figure(fg++); clf; plot(rx_np_log,'+'); title('Octave Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); - figure(fg++); clf; plot(rx_np_log_c,'+'); title('C Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); - - length(rxbuf_in_log_c) - length(rxbuf_in_log) - - stem_sig_and_error(fg++, 111, tx_bits_log_c, tx_bits_log - tx_bits_log_c, 'tx bits', [1 length(tx_bits_log) -1.5 1.5]) - - stem_sig_and_error(fg, 211, real(tx_log_c), real(tx_log - tx_log_c), 'tx re', [1 length(tx_log_c) -0.1 0.1]) - stem_sig_and_error(fg++, 212, imag(tx_log_c), imag(tx_log - tx_log_c), 'tx im', [1 length(tx_log_c) -0.1 0.1]) - - stem_sig_and_error(fg, 211, real(rx_log_c), real(rx_log - rx_log_c), 'rx re', [1 length(rx_log_c) -0.1 0.1]) - stem_sig_and_error(fg++, 212, imag(rx_log_c), imag(rx_log - rx_log_c), 'rx im', [1 length(rx_log_c) -0.1 0.1]) - - %stem_sig_and_error(fg, 211, real(rxbuf_in_log_c), real(rxbuf_in_log - rxbuf_in_log_c), 'rxbuf in re', [1 length(rxbuf_in_log_c) -0.1 0.1]) - %stem_sig_and_error(fg++, 212, imag(rxbuf_in_log_c), imag(rxbuf_in_log - rxbuf_in_log_c), 'rxbuf in im', [1 length(rxbuf_in_log_c) -0.1 0.1]) - - [rxbuf_log_c rxbuf_log] = trim1(rxbuf_log_c,rxbuf_log); - stem_sig_and_error(fg, 211, real(rxbuf_log_c), real(rxbuf_log - rxbuf_log_c), 'rxbuf re', [1 length(rxbuf_log_c) -0.1 0.1]) - stem_sig_and_error(fg++, 212, imag(rxbuf_log_c), imag(rxbuf_log - rxbuf_log_c), 'rxbuf im', [1 length(rxbuf_log_c) -0.1 0.1]) - - [rx_sym_log_c rx_sym_log] = trim2(rx_sym_log_c,rx_sym_log); - stem_sig_and_error(fg, 211, real(rx_sym_log_c), real(rx_sym_log - rx_sym_log_c), 'rx sym re', [1 length(rx_sym_log_c) -1.5 1.5]) - stem_sig_and_error(fg++, 212, imag(rx_sym_log_c), imag(rx_sym_log - rx_sym_log_c), 'rx sym im', [1 length(rx_sym_log_c) -1.5 1.5]) - - % for angles pi and -pi are the same + prx += lnew; + + [rx_bits states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in); + + % log some states for comparison to C + + rxbuf_in_log = [rxbuf_in_log rxbuf_in]; + rxbuf_log = [rxbuf_log states.rxbuf]; + rx_sym_log = [rx_sym_log; states.rx_sym]; + phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log]; + rx_amp_log = [rx_amp_log arx_amp]; + foff_hz_log = [foff_hz_log; states.foff_est_hz]; + timing_est_log = [timing_est_log; states.timing_est]; + sample_point_log = [sample_point_log; states.sample_point]; + rx_np_log = [rx_np_log arx_np]; + rx_bits_log = [rx_bits_log rx_bits]; - [phase_est_pilot_log phase_est_pilot_log_c] = trim2(phase_est_pilot_log, phase_est_pilot_log_c); - d = phase_est_pilot_log - phase_est_pilot_log_c; d = angle(exp(j*d)); - [rx_amp_log rx_amp_log_c] = trim1(rx_amp_log,rx_amp_log_c); - stem_sig_and_error(fg, 211, phase_est_pilot_log_c, d, 'phase est pilot', [1 length(phase_est_pilot_log_c) -1.5 1.5]) - stem_sig_and_error(fg++, 212, rx_amp_log_c, rx_amp_log - rx_amp_log_c, 'rx amp', [1 length(rx_amp_log_c) -1.5 1.5]) - - stem_sig_and_error(fg++, 111, foff_hz_log_c, (foff_hz_log - foff_hz_log_c), 'foff hz', [1 length(foff_hz_log_c) -1.5 1.5]) - - stem_sig_and_error(fg, 211, timing_est_log_c, (timing_est_log - timing_est_log_c), 'timing est', [1 length(timing_est_log_c) -1.5 1.5]) - stem_sig_and_error(fg++, 212, sample_point_log_c, (sample_point_log - sample_point_log_c), 'sample point', [1 length(sample_point_log_c) -1.5 1.5]) - - [rx_bits_log rx_bits_log_c] = trim1(rx_bits_log,rx_bits_log_c); - stem_sig_and_error(fg++, 111, rx_bits_log_c, rx_bits_log - rx_bits_log_c, 'rx bits', [1 length(rx_bits_log) -1.5 1.5]) - - [rxbuf_in_log rxbuf_in_log_c] = trim1(rxbuf_in_log,rxbuf_in_log_c); - [timing_est_log timing_est_log_c] = trim1(timing_est_log,timing_est_log_c'); - [sample_point_log,sample_point_log_c] = trim1(sample_point_log,sample_point_log_c); - [foff_hz_log,foff_hz_log_c] = trim1(foff_hz_log,foff_hz_log_c'); - - % Run through checklist ----------------------------- - pass = true; - pass = check_no_abs(W, W_c, 'W') && pass; - pass = check(tx_bits_log, tx_bits_log_c, 'tx_bits') && pass; - pass = check(tx_log, tx_log_c, 'tx') && pass; - pass = check(rx_log, rx_log_c, 'rx') && pass; - pass = check(rxbuf_in_log, rxbuf_in_log_c, 'rxbuf in') && pass; - pass = check(rxbuf_log, rxbuf_log_c, 'rxbuf') && pass; - pass = check(rx_sym_log, rx_sym_log_c, 'rx_sym') && pass; - pass = check(phase_est_pilot_log, phase_est_pilot_log_c, 'phase_est_pilot', tol=1E-3, its_an_angle=1) && pass; - pass = check(rx_amp_log, rx_amp_log_c, 'rx_amp') && pass; - pass = check(timing_est_log, timing_est_log_c, 'timing_est') && pass; - pass = check(sample_point_log, sample_point_log_c, 'sample_point') && pass; - pass = check(foff_hz_log, foff_hz_log_c, 'foff_est_hz') && pass; - pass = check(rx_bits_log, rx_bits_log_c, 'rx_bits') && pass; - - figure(fg++) - stem(timing_est_log-timing_est_log_c) - figure(fg++) - stem(phase_est_pilot_log-phase_est_pilot_log_c) - end -function [ber berc] = run_ber_test(Nframes,sample_clock_offset_ppm,foff_hz,EbNodB = 100,nheadsamp = 0) - rand("seed",2) - randn("seed",2) - %Nframes = 30; - %sample_clock_offset_ppm = 100; - %foff_hz = 5; - - more off; format; - ofdm_lib; - autotest; - - % --------------------------------------------------------------------- - % Run Octave version - % --------------------------------------------------------------------- - - Ts = 0.018; Tcp = 0.002; Rs = 1/Ts; bps = 2; Nc = 16; Ns = 8; - Rb_real = 1400.0; - states = ofdm_init(bps, Rs, Tcp, Ns, Nc); - ofdm_load_const; - - tx_bits = round(rand(1,Nbitsperframe)); - - % Run tx loop - - tx_bits_log = []; tx_log = []; - for f=1:Nframes - tx_bits_log = [tx_bits_log tx_bits]; - tx_log = [tx_log ofdm_mod(states, tx_bits)]; - end - - % Channel simulation ---------------------------------------------- - - rx_log = sample_clock_offset(tx_log, sample_clock_offset_ppm); - rx_log = freq_shift(rx_log, foff_hz, Fs); - #rx_log(1:nheadsamp) = zeros(nheadsamp,1); - rx_log_l = length(rx_log) - rx_log = [zeros(1,nheadsamp), rx_log]; - rx_log = rx_log(1:rx_log_l); - %rx_log = snr_awgn_channel(rx_log, Fs, Fs, EbNodB); - rx_log = ebno_awgn_channel(rx_log, Fs, Rb_real, EbNodB); +% --------------------------------------------------------------------- +% Run C version and plot Octave and C states and differences +% --------------------------------------------------------------------- - rx_vec = fopen("tofdm_rx_vec","wb+"); - fwrite(rx_vec,im_re_interleave(rx_log),"single"); - fclose(rx_vec); +% Override default path by setting path_to_tofdm = "/your/path/to/tofdm" - % Rx --------------------------------------------------------------- +if exist("path_to_tofdm", "var") == 0 + path_to_tofdm = "../build_linux/unittest/tofdm"; +end +system(path_to_tofdm); - % Init rx with ideal timing so we can test with timing estimation disabled +load tofdm_out.txt; - Nsam = length(rx_log); - prx = 1; - nin = Nsamperframe+2*(M+Ncp); - %states.rxbuf(Nrxbuf-nin+1:Nrxbuf) = rx_log(prx:nin); - %prx += nin; +fg = 1; +figure(fg++); clf; plot(rx_np_log,'+'); title('Octave Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); +figure(fg++); clf; plot(rx_np_log_c,'+'); title('C Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); - rxbuf_log = []; rxbuf_in_log = []; rx_sym_log = []; foff_hz_log = []; - timing_est_log = []; sample_point_log = []; - phase_est_pilot_log = []; rx_amp_log = []; - rx_np_log = []; rx_bits_log = []; +stem_sig_and_error(fg++, 111, tx_bits_log_c, tx_bits_log - tx_bits_log_c, 'tx bits', [1 length(tx_bits_log) -1.5 1.5]) - states.timing_en = 1; - states.foff_est_en = 1; - states.phase_est_en = 1; +stem_sig_and_error(fg, 211, real(tx_log_c), real(tx_log - tx_log_c), 'tx re', [1 length(tx_log_c) -0.1 0.1]) +stem_sig_and_error(fg++, 212, imag(tx_log_c), imag(tx_log - tx_log_c), 'tx im', [1 length(tx_log_c) -0.1 0.1]) - if states.timing_en == 0 - % manually set ideal timing instant - states.sample_point = Ncp; - end +stem_sig_and_error(fg, 211, real(rx_log_c), real(rx_log - rx_log_c), 'rx re', [1 length(rx_log_c) -0.1 0.1]) +stem_sig_and_error(fg++, 212, imag(rx_log_c), imag(rx_log - rx_log_c), 'rx im', [1 length(rx_log_c) -0.1 0.1]) - have_frame = true; - while have_frame - % insert samples at end of buffer, set to zero if no samples - % available to disable phase estimation on future pilots on last - % frame of simulation - - nin = states.nin; - lnew = min(Nsam-prx+1,nin); - rxbuf_in = zeros(1,nin); - %printf("nin: %d prx: %d lnew: %d\n", nin, prx, lnew); - if lnew - rxbuf_in(1:lnew) = rx_log(prx:prx+lnew-1); - end - prx += nin; - if prx > Nsam - have_frame = false; - end - - [rx_bits states aphase_est_pilot_log arx_np arx_amp] = ofdm_demod(states, rxbuf_in); - % log some states for comparison to C - - rxbuf_in_log = [rxbuf_in_log rxbuf_in]; - rxbuf_log = [rxbuf_log states.rxbuf]; - rx_sym_log = [rx_sym_log; states.rx_sym]; - phase_est_pilot_log = [phase_est_pilot_log; aphase_est_pilot_log]; - rx_amp_log = [rx_amp_log arx_amp]; - foff_hz_log = [foff_hz_log; states.foff_est_hz]; - timing_est_log = [timing_est_log; states.timing_est]; - sample_point_log = [sample_point_log; states.sample_point]; - rx_np_log = [rx_np_log arx_np]; - rx_bits_log = [rx_bits_log rx_bits]; - - end +stem_sig_and_error(fg, 211, real(rxbuf_in_log_c), real(rxbuf_in_log - rxbuf_in_log_c), 'rxbuf in re', [1 length(rxbuf_in_log_c) -0.1 0.1]) +stem_sig_and_error(fg++, 212, imag(rxbuf_in_log_c), imag(rxbuf_in_log - rxbuf_in_log_c), 'rxbuf in im', [1 length(rxbuf_in_log_c) -0.1 0.1]) - % --------------------------------------------------------------------- - % Run C version and plot Octave and C states and differences - % --------------------------------------------------------------------- +stem_sig_and_error(fg, 211, real(rxbuf_log_c), real(rxbuf_log - rxbuf_log_c), 'rxbuf re', [1 length(rxbuf_log_c) -0.1 0.1]) +stem_sig_and_error(fg++, 212, imag(rxbuf_log_c), imag(rxbuf_log - rxbuf_log_c), 'rxbuf im', [1 length(rxbuf_log_c) -0.1 0.1]) - % Override default path by setting path_to_tofdm = "/your/path/to/tofdm" +stem_sig_and_error(fg, 211, real(rx_sym_log_c), real(rx_sym_log - rx_sym_log_c), 'rx sym re', [1 length(rx_sym_log_c) -1.5 1.5]) +stem_sig_and_error(fg++, 212, imag(rx_sym_log_c), imag(rx_sym_log - rx_sym_log_c), 'rx sym im', [1 length(rx_sym_log_c) -1.5 1.5]) - if exist("path_to_tofdm", "var") == 0 - path_to_tofdm = "../build_linux/unittest/tofdm"; - end - system(path_to_tofdm); +% for angles pi and -pi are the same - load tofdm_out.txt; - % Generated with modem probe thing - load ofdm_test.txt; +d = phase_est_pilot_log - phase_est_pilot_log_c; d = angle(exp(j*d)); - front_cutoff = 10000; - rx_np_log = rx_np_log(front_cutoff/2:length(rx_np_log)); - rx_np_log_c = rx_np_log_c(front_cutoff/2:length(rx_np_log_c)); +stem_sig_and_error(fg, 211, phase_est_pilot_log_c, d, 'phase est pilot', [1 length(phase_est_pilot_log_c) -1.5 1.5]) +stem_sig_and_error(fg++, 212, rx_amp_log_c, rx_amp_log - rx_amp_log_c, 'rx amp', [1 length(rx_amp_log_c) -1.5 1.5]) - fg = 1; - figure(fg++); clf; plot(foff_hz_log_c); title('foff hz C'); - figure(fg++); clf; plot(rx_np_log,'+'); title('Octave Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); - figure(fg++); clf; plot(rx_np_log_c,'+'); title('C Scatter Diagram'); axis([-1.5 1.5 -1.5 1.5]); - +stem_sig_and_error(fg++, 111, foff_hz_log_c, (foff_hz_log - foff_hz_log_c), 'foff hz', [1 length(foff_hz_log_c) -1.5 1.5]) - [rx_bits_log rx_bits_log_c] = trim1(rx_bits_log,rx_bits_log_c); - [tx_bits_log rx_bits_log] = trim1(tx_bits_log,rx_bits_log); - [tx_bits_log rx_bits_log_c] = trim1(tx_bits_log,rx_bits_log_c); - - bitcnt = length(rx_bits_log_c); - rx_bits = rx_bits_log_c(front_cutoff:length(rx_bits_log_c)); - tx_bits = tx_bits_log(front_cutoff:length(tx_bits_log)); - - ber = 1; - for offset = (1:100) - nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset))); - bern = nerr/(bitcnt-offset); - if(bern < ber) - ox = offset; - best_nerr = nerr; - end - ber = min([ber bern]); - end - offset = ox; - berc = ber; - - figure(fg++) - plot(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset))) - title('C bit errors') - - rx_bits = rx_bits_log(front_cutoff:length(rx_bits)); - ber = 1; - - for offset = (1:100) - nerr = sum(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset))); - bern = nerr/(bitcnt-offset); - if(bern < ber) - ox = offset; - best_nerr = nerr; - end - ber = min([ber bern]); - end - offset = ox; - ber = ber; +stem_sig_and_error(fg, 211, timing_est_log_c, (timing_est_log - timing_est_log_c), 'timing est', [1 length(timing_est_log_c) -1.5 1.5]) +stem_sig_and_error(fg++, 212, sample_point_log_c, (sample_point_log - sample_point_log_c), 'sample point', [1 length(sample_point_log_c) -1.5 1.5]) - figure(fg++) - plot(xor(rx_bits(offset:length(rx_bits)),tx_bits(1:length(rx_bits)+1-offset))) - title('Octave bit errors') - - - figure(fg++) - plot(20*log10(abs(fft(rx_log))), (1:length(rx_log)) .* (Fs/length(rx_log) )) +stem_sig_and_error(fg++, 111, rx_bits_log_c, rx_bits_log - rx_bits_log_c, 'rx bits', [1 length(rx_bits_log) -1.5 1.5]) -endfunction +% Run through checklist ----------------------------- -% FLTK plots work better than gnuplot or the default on my machine -- Brady -try - graphics_toolkit('fltk') -catch err - warning(err.identifier, err.message); -end_try_catch +check_no_abs(W, W_c, 'W'); +check(tx_bits_log, tx_bits_log_c, 'tx_bits'); +check(tx_log, tx_log_c, 'tx'); +check(rx_log, rx_log_c, 'rx'); +check(rxbuf_in_log, rxbuf_in_log_c, 'rxbuf in'); +check(rxbuf_log, rxbuf_log_c, 'rxbuf'); +check(rx_sym_log, rx_sym_log_c, 'rx_sym'); +check(phase_est_pilot_log, phase_est_pilot_log_c, 'phase_est_pilot', tol=1E-3, its_an_angle=1); +check(rx_amp_log, rx_amp_log_c, 'rx_amp'); +check(timing_est_log, timing_est_log_c, 'timing_est'); +check(sample_point_log, sample_point_log_c, 'sample_point'); +check(foff_hz_log, foff_hz_log_c, 'foff_est_hz'); +check(rx_bits_log, rx_bits_log_c, 'rx_bits'); -%run_ofdm_test(60,100,60,5,10000) diff --git a/codec2-dev/src/codec2_ofdm.h b/codec2-dev/src/codec2_ofdm.h index fc1833bd..b335df88 100644 --- a/codec2-dev/src/codec2_ofdm.h +++ b/codec2-dev/src/codec2_ofdm.h @@ -55,8 +55,8 @@ void ofdm_destroy(struct OFDM *); void ofdm_mod(struct OFDM *, COMP *, const int *); void ofdm_demod(struct OFDM *, int *, COMP *); int ofdm_get_nin(struct OFDM *); -int ofdm_get_samples_per_frame(struct OFDM *ofdm); -int ofdm_get_max_samples_per_frame(struct OFDM *ofdm); +int ofdm_get_samples_per_frame(void); +int ofdm_get_max_samples_per_frame(void); /* option setters */ diff --git a/codec2-dev/src/ofdm.c b/codec2-dev/src/ofdm.c index 9229b176..52079af4 100644 --- a/codec2-dev/src/ofdm.c +++ b/codec2-dev/src/ofdm.c @@ -38,39 +38,10 @@ #include "ofdm_internal.h" #include "codec2_ofdm.h" #include "kiss_fft.h" -#include "modem_probe.h" - -/* - * Optimization TODOs - * CPU: - * [ ] Remove cexpf() from loops (can be replaced with single cexpf() and multiplication in loop) - * [ ] Look for/fix potential aliasing slowdowns (ie using array[i] as accumulator) - * [ ] Possibly dd a bunch of const arrays and maybe figure out restrict - * [ ] Run a bit of profiling on ofdm - * - * Mem: - * [ ] Make RX debug logs optional or modem probes - */ /* Concrete definition of 700D parameters */ -const struct OFDM_CONFIG OFDM_CONFIG_700D_C = { - .Nc = 16, - .Ts = 0.018f, - .Rs = (1.0/0.018), - .Fs = 8000, - .bps = 2, - .Tcp = -1, - .Ns = 8, - .M = 144, - .Ncp = 16, - .Fcenter = 1500, - .FtWindowWidth = 11, - .BitsPerFrame = 224, - .SampsPerFrame = 1280, - .SampsPerFrameMax = 1320, - .RxBufSize = 4320, - .RowsPerFrame = 7 -}; +const struct OFDM_CONFIG OFDM_CONFIG_700D_C = +{.a = 0}; /* Pointer to 700D config */ const struct OFDM_CONFIG * OFDM_CONFIG_700D = &OFDM_CONFIG_700D_C; @@ -82,9 +53,8 @@ static void idft(struct OFDM *, complex float *, complex float *); static complex float vector_sum(complex float *, int); static complex float qpsk_mod(int *); static void qpsk_demod(complex float, int *); -static void ofdm_txframe(struct OFDM *, complex float [], complex float *); -//static int coarse_sync(struct OFDM *, complex float *, int); -static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length, float *foff_out, float *ratio_out); +static void ofdm_txframe(struct OFDM *, complex float [OFDM_SAMPLESPERFRAME], complex float *); +static int coarse_sync(struct OFDM *, complex float *, int); /* Defines */ @@ -136,14 +106,12 @@ static void qpsk_demod(complex float symbol, int *bits) { static void idft(struct OFDM *ofdm, complex float *result, complex float *vector) { int row, col; - int Nc = ofdm->config.Nc; - int M = ofdm->config.M; - for (row = 0; row < M; row++) { + for (row = 0; row < OFDM_M; row++) { result[row] = 0.0f + 0.0f * I; - for (col = 0; col < (Nc + 2); col++) { - result[row] = result[row] + (vector[col] * (ofdm->W[col][row] / (float) M)); /* complex result */ + for (col = 0; col < (OFDM_NC + 2); col++) { + result[row] = result[row] + (vector[col] * (ofdm->W[col][row] / (float) OFDM_M)); /* complex result */ } } } @@ -152,13 +120,11 @@ static void idft(struct OFDM *ofdm, complex float *result, complex float *vector static void dft(struct OFDM *ofdm, complex float *result, complex float *vector) { int row, col; - int Nc = ofdm->config.Nc; - int M = ofdm->config.M; - for (col = 0; col < (Nc + 2); col++) { + for (col = 0; col < (OFDM_NC + 2); col++) { result[col] = 0.0f + 0.0f * I; - for (row = 0; row < M; row++) { + for (row = 0; row < OFDM_M; row++) { result[col] = result[col] + (vector[row] * conjf(ofdm->W[col][row])); /* complex result */ } } @@ -176,225 +142,52 @@ static complex float vector_sum(complex float *a, int num_elements) { return sum; } -static int coarser_sync(struct OFDM *ofdm, complex float *rx, int length, float *foff_out, float *ratio_out) { - complex float csam; - const int Fs = ofdm->config.Fs; - const int M = ofdm->config.M; - const int Ncp = ofdm->config.Ncp; - const int SampsPerFrame = ofdm->config.SampsPerFrame; - const int Ncorr = length - (M + Ncp); - - float corr[Ncorr]; - float corr_t; - int NPSamp = M + Ncp; - int i, j; - float max_corr = 0; - float max_mag = 0; - float mag = 0.0f; - int t_est = 0; - int fmax = 20; - int pmax_i,nmax_i; - float pmax,nmax; - float foff_est; - - for (i = 0; i < Ncorr; i++) { - complex float temp = 0.0f + 0.0f * I; - - for (j = 0; j < (M + Ncp); j++) { - csam = conjf(ofdm->pilot_samples[j]); - temp = temp + (rx[i + j] * csam); - //temp = temp + (rx[i + j + SampsPerFrame] * csam); - } - - corr_t = cabsf(temp); - /* Magnitude of incoming samples */ - //mag = max( cabsf(rx[i]) , cabsf(rx[i+SampsPerFrame]) ); - mag = cabsf(rx[i]); - max_mag = max(mag, max_mag); - - /* find the max magnitude and its index */ - if(corr_t > max_corr){ - max_corr = corr_t; - t_est = i; - } - } - - if(ratio_out != NULL) - *ratio_out = max_corr/max_mag; - //fprintf(stderr,"Ratio is %f; mag is %f corr is %f\n",max_corr/max_mag,max_mag,max_corr); - - /* Skip coarse frequency est if we don't have anywhere to put result */ - if(foff_out == NULL){ - return t_est; - } - - /* Coarse frequency estimation */ - kiss_fft_cfg fftcfg = ofdm->sync_fft_cfg; - complex float fft_in[Fs]; - complex float fft_out[Fs]; - float C[Fs]; - /* Zero FFT input array */ - for(i = 0; i < Fs; i++){ - fft_in[i] = 0; - } - /* shift and copy in NPsam samples to front of buffer for FFT'ing */ - for(i = 0; i < NPSamp; i++){ - fft_in[i] = rx[i + t_est] * conjf(ofdm->pilot_samples[i]); - } - kiss_fft(fftcfg,(kiss_fft_cpx*)fft_in,(kiss_fft_cpx*)fft_out); - - /* Copy into output corr array, taking magnitude */ - /* TODO: May be able to skip sqrt() in abs, since we're just looking for a maximum point down the line */ - for(i = 0; i < Fs; i++){ - C[i] = cabsf(fft_out[i]); - } - - // /* shift and copy in NPsam samples to front of buffer for FFT'ing */ - // for(i = 0; i < NPSamp; i++){ - // fft_in[i] = rx[i + t_est + SampsPerFrame] * conjf(ofdm->pilot_samples[i]); - // } - // kiss_fft(fftcfg,(kiss_fft_cpx*)fft_in,(kiss_fft_cpx*)fft_out); - - // /* Add into output corr array */ - // for(i = 0; i < Fs; i++){ - // C[i] += cabsf(fft_out[i]); - // } - - pmax_i = 0; - nmax_i = Fs; - pmax = nmax = 0; - /* Search the positive and negative sides of the FFT to +/-fmax for peak */ - for(i = 0; i < fmax; i++){ - if(C[i] > pmax){ - pmax = C[i]; - pmax_i = i; - } - } - for(i = Fs-1; i > Fs-fmax; i--){ - if(C[i] > nmax){ - nmax = C[i]; - nmax_i = i; - } - } - foff_est = (pmax > nmax) ? pmax_i : (nmax_i-Fs+1); - - //fprintf(stderr,"foff_est is %f, tpos is %d, ratio is %f\n",foff_est,t_est,max_corr/max_mag); - - *foff_out = foff_est; - - return t_est; -} - /* * Correlates the OFDM pilot symbol samples with a window of received * samples to determine the most likely timing offset. Combines two * frames pilots so we need at least Nsamperframe+M+Ncp samples in rx. */ -static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length, float *foff_out, float *ratio_out) { +static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length) { complex float csam; - const int Fs = ofdm->config.Fs; - const int M = ofdm->config.M; - const int Ncp = ofdm->config.Ncp; - const int SampsPerFrame = ofdm->config.SampsPerFrame; - const int Ncorr = length - (SampsPerFrame + (M + Ncp)); - + int Ncorr = length - (OFDM_SAMPLESPERFRAME + (OFDM_M + OFDM_NCP)); + int Fs = OFDM_FS; + int SFrame = OFDM_SAMPLESPERFRAME; float corr[Ncorr]; - float corr_t; - int NPSamp = M + Ncp; int i, j; - float max_corr = 0; - float max_mag = 0; - float mag = 0.0f; - int t_est = 0; - int fmax = 20; - int pmax_i,nmax_i; - float pmax,nmax; - float foff_est; for (i = 0; i < Ncorr; i++) { complex float temp = 0.0f + 0.0f * I; - for (j = 0; j < (M + Ncp); j++) { + for (j = 0; j < (OFDM_M + OFDM_NCP); j++) { csam = conjf(ofdm->pilot_samples[j]); temp = temp + (rx[i + j] * csam); - temp = temp + (rx[i + j + SampsPerFrame] * csam); + temp = temp + (rx[i + j + SFrame] * csam); } - corr_t = cabsf(temp); - /* Magnitude of incoming samples */ - mag = max( cabsf(rx[i]) , cabsf(rx[i+SampsPerFrame]) ); - max_mag = max(mag, max_mag); - - /* find the max magnitude and its index */ - if(corr_t > max_corr){ - max_corr = corr_t; - t_est = i; - } + corr[i] = cabsf(temp); } - if(ratio_out != NULL){ - *ratio_out = (max_corr + 1e-9)/(max_mag + 1e-9); - } + /* find the max magnitude and its index */ + + float mag = 0.0f; + int t_est = 0; - /* Skip coarse frequency est if we don't have anywhere to put result */ - if(foff_out == NULL){ - return t_est; + for (i = 0; i < Ncorr; i++) { + if (corr[i] > mag) { + mag = corr[i]; + t_est = i; + } } /* Coarse frequency estimation */ - kiss_fft_cfg fftcfg = ofdm->sync_fft_cfg; + /* TODO: Move FFT config to ofdm init and ofdm struct */ + kiss_fft_cfg fftcfg = kiss_fft_alloc(Fs,0,NULL,NULL); complex float fft_in[Fs]; complex float fft_out[Fs]; - float C[Fs]; - /* Zero FFT input array */ - for(i = 0; i < Fs; i++){ - fft_in[i] = 0; - } - /* shift and copy in NPsam samples to front of buffer for FFT'ing */ - for(i = 0; i < NPSamp; i++){ - fft_in[i] = rx[i + t_est] * conjf(ofdm->pilot_samples[i]); - } - kiss_fft(fftcfg,(kiss_fft_cpx*)fft_in,(kiss_fft_cpx*)fft_out); - - /* Copy into output corr array, taking magnitude */ - /* TODO: May be able to skip sqrt() in abs, since we're just looking for a maximum point down the line */ - for(i = 0; i < Fs; i++){ - C[i] = cabsf(fft_out[i]); - } - - /* shift and copy in NPsam samples to front of buffer for FFT'ing */ - for(i = 0; i < NPSamp; i++){ - fft_in[i] = rx[i + t_est + SampsPerFrame] * conjf(ofdm->pilot_samples[i]); - } - kiss_fft(fftcfg,(kiss_fft_cpx*)fft_in,(kiss_fft_cpx*)fft_out); - - /* Add into output corr array */ - for(i = 0; i < Fs; i++){ - C[i] += cabsf(fft_out[i]); - } - - pmax_i = 0; - nmax_i = Fs; - pmax = nmax = 0; - /* Search the positive and negative sides of the FFT to +/-fmax for peak */ - for(i = 0; i < fmax; i++){ - if(C[i] > pmax){ - pmax = C[i]; - pmax_i = i; - } - } - for(i = Fs-1; i > Fs-fmax; i--){ - if(C[i] > nmax){ - nmax = C[i]; - nmax_i = i; - } - } - foff_est = (pmax > nmax) ? pmax_i : (nmax_i-Fs+1); - //foff_est = pmax_i; - //fprintf(stderr,"foff_est is %f, tpos is %d, ratio is %f\n",foff_est,t_est%1280,max_corr/max_mag); + + - *foff_out = foff_est; return t_est; } @@ -404,119 +197,64 @@ static int coarse_sync(struct OFDM *ofdm, complex float *rx, int length, float * * ---------------------------------------------- */ -static void ofdm_txframe(struct OFDM *ofdm, complex float tx[], +static void ofdm_txframe(struct OFDM *ofdm, complex float tx[OFDM_SAMPLESPERFRAME], complex float *tx_sym_lin) { - - const int Nc = ofdm->config.Nc; - const int M = ofdm->config.M; - const int Ncp = ofdm->config.Ncp; - const int Ns = ofdm->config.Ns; - const int RowsPerFrame = ofdm->config.RowsPerFrame; - - complex float aframe[Ns][Nc + 2]; - complex float asymbol[M]; - complex float asymbol_cp[M + Ncp]; + complex float aframe[OFDM_NS][OFDM_NC + 2]; + complex float asymbol[OFDM_M]; + complex float asymbol_cp[OFDM_M + OFDM_NCP]; int i, j, k, m; /* initialize aframe to complex zero */ - for (i = 0; i < Ns; i++) { - for (j = 0; j < (Nc + 2); j++) { + for (i = 0; i < OFDM_NS; i++) { + for (j = 0; j < (OFDM_NC + 2); j++) { aframe[i][j] = 0.0f + 0.0f * I; } } /* copy in a row of complex pilots to first row */ - for (i = 0; i < (Nc + 2); i++) { + for (i = 0; i < (OFDM_NC + 2); i++) { aframe[0][i] = ofdm->pilots[i]; } /* Place symbols in multi-carrier frame with pilots */ /* This will place boundary values of complex zero around data */ - for (i = 1; i <= RowsPerFrame; i++) { + for (i = 1; i <= OFDM_ROWSPERFRAME; i++) { /* copy in the Nc complex values with [0 Nc 0] or (Nc + 2) total */ - for (j = 1; j < (Nc + 1); j++) { - aframe[i][j] = tx_sym_lin[((i - 1) * Nc) + (j - 1)]; + for (j = 1; j < (OFDM_NC + 1); j++) { + aframe[i][j] = tx_sym_lin[((i - 1) * OFDM_NC) + (j - 1)]; } } /* OFDM up-convert symbol by symbol so we can add CP */ - for (i = 0, m = 0; i < Ns; i++, m += (M + Ncp)) { + for (i = 0, m = 0; i < OFDM_NS; i++, m += (OFDM_M + OFDM_NCP)) { idft(ofdm, asymbol, aframe[i]); /* Copy the last Ncp columns to the front */ - for (j = (M - Ncp), k = 0; j < M; j++, k++) { + for (j = (OFDM_M - OFDM_NCP), k = 0; j < OFDM_M; j++, k++) { asymbol_cp[k] = asymbol[j]; } /* Now copy the whole row after it */ - for (j = Ncp, k = 0; k < M; j++, k++) { + for (j = OFDM_NCP, k = 0; k < OFDM_M; j++, k++) { asymbol_cp[j] = asymbol[k]; } /* Now move row to the tx reference */ - for (j = 0; j < (M + Ncp); j++) { + for (j = 0; j < (OFDM_M + OFDM_NCP); j++) { tx[m + j] = asymbol_cp[j]; } } } -/* - * Utility function to allocate a 2d array of dimension - * [sx][sy] with elements of elem size - */ -void **alloc_doubleary(size_t sx,size_t sy,size_t elem){ - size_t i; - /* Allocate pointer array */ - char ** ary = malloc(sizeof(void*) * sx); - if(ary == NULL){ - return NULL; - } - /* Allocate data array. */ - *ary = malloc(elem * sx * sy); - if(*ary == NULL){ - free(ary); - return NULL; - } - /* Fill out pointer array to point to rows of data array */ - for(i=0; iconfig,(void*)config,sizeof(struct OFDM_CONFIG)); - const int Fs = ofdm->config.Fs; - const float Rs = ofdm->config.Rs; - const int Nc = ofdm->config.Nc; - const int M = ofdm->config.M; - const int Ncp = ofdm->config.Ncp; - const int Ns = ofdm->config.Ns; - const int RowsPerFrame = ofdm->config.RowsPerFrame; - const int SampsPerFrame = ofdm->config.SampsPerFrame; - const int Fcenter = ofdm->config.Fcenter; - complex float temp[M]; - - /* Allocate various buffers */ - /* TODO: cleanup after failed malloc */ - ofdm->pilot_samples = malloc(sizeof(complex float) * (M + Ncp)); - if(ofdm->pilot_samples == NULL) goto ofdm_fail_cleanup; - - ofdm->pilots = malloc(sizeof(complex float) * (Nc+2)); - if(ofdm->pilots == NULL) goto ofdm_fail_cleanup; - - ofdm->rxbuf = malloc(sizeof(complex float) * ofdm->config.RxBufSize); - if(ofdm->rxbuf == NULL) goto ofdm_fail_cleanup; - - ofdm->coarse_rxbuf = malloc(sizeof(complex float) * ofdm->config.RxBufSize); - if(ofdm->coarse_rxbuf == NULL) goto ofdm_fail_cleanup; - - ofdm->w = malloc(sizeof(float) * (Nc+2)); - if(ofdm->w == NULL) goto ofdm_fail_cleanup; - - ofdm->rx_amp = malloc(sizeof(float) * RowsPerFrame * Nc); - if(ofdm->rx_amp == NULL) goto ofdm_fail_cleanup; - - ofdm->aphase_est_pilot_log = malloc(sizeof(float) * RowsPerFrame * Nc); - if(ofdm->aphase_est_pilot_log == NULL) goto ofdm_fail_cleanup; - - ofdm->sync_fft_cfg = kiss_fft_alloc(Fs,0,NULL,NULL); - if(ofdm->sync_fft_cfg == NULL) goto ofdm_fail_cleanup; /* store complex BPSK pilot symbols */ - for (i = 0; i < (Nc + 2); i++) { + for (i = 0; i < (OFDM_NC + 2); i++) { ofdm->pilots[i] = ((float) pilotvalues[i]) + 0.0f * I; } /* carrier tables for up and down conversion */ - float Fcf = (float)Fcenter; - float Ncf = (float)Nc; - float Fsf = (float)Fs; - int Nlower = floorf(( Fcf - Rs * (Ncf / 2)) / Rs); - for (i = 0, j = Nlower; i < (Nc + 2); i++, j++) { + int Nlower = floorf((OFDM_CENTRE - OFDM_RS * (OFDM_NC / 2)) / OFDM_RS); + + for (i = 0, j = Nlower; i < (OFDM_NC + 2); i++, j++) { /* * 2 * pi * j/144 j=19..36 * j = 1 kHz to 2 kHz (1.5 kHz center) */ - ofdm->w[i] = TAU * (float) j / (Fsf / Rs); + ofdm->w[i] = TAU * (float) j / (OFDM_FS / OFDM_RS); } - ofdm->W = (complex float **) alloc_doubleary(Nc+2,M,sizeof(complex float)); - if(ofdm->W == NULL) goto ofdm_fail_cleanup; - - ofdm->rx_sym = (complex float **) alloc_doubleary(Ns+3,Nc+2,sizeof(complex float)); - if(ofdm->rx_sym == NULL) goto ofdm_fail_cleanup; - - for (i = 0; i < (Nc + 2); i++) { - - for (j = 0; j < M; j++) { + for (i = 0; i < (OFDM_NC + 2); i++) { + for (j = 0; j < OFDM_M; j++) { ofdm->W[i][j] = cexpf(I * ofdm->w[i] * j); } } - for (i = 0; i < (Ns + 3); i++) { - for (j = 0; j < (Nc + 2); j++) { + for (i = 0; i < (OFDM_NS + 3); i++) { + for (j = 0; j < (OFDM_NC + 2); j++) { ofdm->rx_sym[i][j] = 0.0f + 0.0f * I; } } - - for(i = 0; i < ofdm->config.RxBufSize; i++){ - ofdm->rxbuf[i] = 0; - } /* default settings of options and states */ @@ -628,12 +317,12 @@ struct OFDM *ofdm_create(const struct OFDM_CONFIG * config) { ofdm->foff_est_hz = 0.0f; ofdm->sample_point = 0; ofdm->timing_est = 0; - ofdm->nin = SampsPerFrame; - ofdm->sync_count = 0; - ofdm->frame_point = 0; + ofdm->nin = OFDM_SAMPLESPERFRAME; /* create the OFDM waveform */ + complex float temp[OFDM_M]; + idft(ofdm, temp, ofdm->pilots); /* @@ -642,65 +331,33 @@ struct OFDM *ofdm_create(const struct OFDM_CONFIG * config) { /* first copy the last Cyclic Prefix (CP) values */ - for (i = 0, j = (M - Ncp); i < Ncp; i++, j++) { + for (i = 0, j = (OFDM_M - OFDM_NCP); i < OFDM_NCP; i++, j++) { ofdm->pilot_samples[i] = temp[j]; } /* Now copy the whole thing after the above */ - for (i = Ncp, j = 0; j < M; i++, j++) { + for (i = OFDM_NCP, j = 0; j < OFDM_M; i++, j++) { ofdm->pilot_samples[i] = temp[j]; } return ofdm; /* Success */ - - /* Deal with a failed malloc or something */ - ofdm_fail_cleanup: - /* Make sure ofdm points to a vaild struct */ - if(ofdm == NULL) - return NULL; - - free_nc(ofdm->pilots); - free_nc(ofdm->pilot_samples); - free_nc(ofdm->rxbuf); - free_nc(ofdm->coarse_rxbuf); - free_nc(ofdm->w); - free_nc(ofdm->rx_amp); - free_nc(ofdm->aphase_est_pilot_log); - free_doubleary((void**)ofdm->W); - free_doubleary((void**)ofdm->rx_sym); - if(ofdm->sync_fft_cfg!=NULL) kiss_fft_free(ofdm->sync_fft_cfg); - free_nc(ofdm); - - return NULL; } void ofdm_destroy(struct OFDM *ofdm) { - if(ofdm == NULL) - return; - - free_nc(ofdm->pilots); - free_nc(ofdm->pilot_samples); - free_nc(ofdm->rxbuf); - free_nc(ofdm->w); - free_nc(ofdm->rx_amp); - free_nc(ofdm->aphase_est_pilot_log); - free_doubleary((void**)ofdm->W); - free_doubleary((void**)ofdm->rx_sym); - if(ofdm->sync_fft_cfg!=NULL) kiss_fft_free(ofdm->sync_fft_cfg); - free_nc(ofdm); + free(ofdm); } int ofdm_get_nin(struct OFDM *ofdm) { return ofdm->nin; } -int ofdm_get_samples_per_frame(struct OFDM *ofdm) { - return ofdm->config.SampsPerFrame; +int ofdm_get_samples_per_frame() { + return OFDM_SAMPLESPERFRAME; } -int ofdm_get_max_samples_per_frame(struct OFDM *ofdm) { - return ofdm->config.SampsPerFrameMax; +int ofdm_get_max_samples_per_frame() { + return OFDM_MAX_SAMPLESPERFRAME; } void ofdm_set_verbose(struct OFDM *ofdm, int level) { @@ -712,7 +369,7 @@ void ofdm_set_timing_enable(struct OFDM *ofdm, bool val) { if (ofdm->timing_en == false) { /* manually set ideal timing instant */ - ofdm->sample_point = (ofdm->config.Ncp - 1); + ofdm->sample_point = (OFDM_NCP - 1); } } @@ -734,30 +391,20 @@ void ofdm_set_off_est_hz(struct OFDM *ofdm, float val) { * -------------------------------------- */ -void ofdm_mod(struct OFDM *ofdm, COMP result[], const int *tx_bits) { - - const int Bps = ofdm->config.bps; - const int SampsPerFrame = ofdm->config.SampsPerFrame; - const int BitsPerFrame = ofdm->config.BitsPerFrame; - - int length = BitsPerFrame / Bps; - - - /* Note: COMP and complex float have the same memory layout */ - /* see iso/iec 9899:1999 sec 6.2.5.13 */ - complex float * tx = (complex float *)&result[0]; - +void ofdm_mod(struct OFDM *ofdm, COMP result[OFDM_SAMPLESPERFRAME], const int *tx_bits) { + int length = OFDM_BITSPERFRAME / OFDM_BPS; + complex float tx[OFDM_SAMPLESPERFRAME]; complex float tx_sym_lin[length]; int dibit[2]; int s, i; - if (Bps == 1) { + if (OFDM_BPS == 1) { /* Here we will have Nbitsperframe / 1 */ for (s = 0; s < length; s++) { tx_sym_lin[s] = (float) (2 * tx_bits[s] - 1) + 0.0f * I; } - } else if (Bps == 2) { + } else if (OFDM_BPS == 2) { /* Here we will have Nbitsperframe / 2 */ for (s = 0, i = 0; i < length; s += 2, i++) { @@ -770,139 +417,13 @@ void ofdm_mod(struct OFDM *ofdm, COMP result[], const int *tx_bits) { ofdm_txframe(ofdm, tx, tx_sym_lin); /* convert to comp */ - /* Note: COMP and complex float have the same memory layout */ - /* see iso/iec 9899:1999 sec 6.2.5.13 */ - /*for (i = 0; i < SampsPerFrame; i++) { + + for (i = 0; i < OFDM_SAMPLESPERFRAME; i++) { result[i].real = crealf(tx[i]); result[i].imag = cimagf(tx[i]); - }*/ -} - -/* - * Like ofdm_demod, but handles sync and large frame offset - */ -//static int coarser_sync(struct OFDM *ofdm, complex float *rx, int length, float *foff_out, float *ratio_out) -void ofdm_demod_coarse(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in){ - const int SampsPerFrame = ofdm->config.SampsPerFrame; - const int RxBufSize = ofdm->config.RxBufSize; - const int Fs = ofdm->config.Fs; - const int M = ofdm->config.M; - const int Ncp = ofdm->config.Ncp; - - int sync = ofdm->sync_count; - int frame_point = ofdm->frame_point; - int nin = ofdm->nin; - complex float * coarse_rxbuf = ofdm->coarse_rxbuf; - complex float rxbuf_temp[RxBufSize]; - - - /* Copy in nin samples to end of buffer */ - memmove(&coarse_rxbuf[0],&coarse_rxbuf[nin], (RxBufSize - nin) * sizeof(coarse_rxbuf[0])); - memcpy (&coarse_rxbuf[RxBufSize - nin],rxbuf_in, nin * sizeof(coarse_rxbuf[0])); - - ofdm_demod(ofdm,rx_bits,&coarse_rxbuf[RxBufSize - 2*SampsPerFrame + frame_point]); - - int frame_pos_est; - float fshift; - float corr_ratio_max = 0; - float corr_ratio; - float foff_out; - float fshift_max; - int frame_pos_max = 0; - /* NCO for freq. shifting */ - complex float shift_nco = 1; - complex float shift_nco_dph = 0; - - if(sync <= 0 ) { - /* Do coarse search over frequency range in 20hz slices*/ - for(fshift = -103; fshift < 100; fshift+=40){ - //fshift = 0;{ - /* Shift block of samples by fshift */ - shift_nco_dph = cexpf(2*M_PI*(fshift/(float)Fs)*I); - shift_nco = 1; - for(int i=0; i corr_ratio_max){ - corr_ratio_max = corr_ratio; - fshift_max = fshift; - frame_pos_max = frame_pos_est; - } - } - frame_pos_est = (frame_pos_max % SampsPerFrame) - (M+Ncp) + 1; - if(corr_ratio_max > 0.5 && sync <=0){ - /* Now that we think we can sync, do a freq. domain search in the smaller range for our signal */ - shift_nco_dph = cexpf(2*M_PI*(fshift_max/(float)Fs)*I); - shift_nco = 1; - for(int i=0; isync_count = 2; - ofdm->frame_point = frame_pos_est; - ofdm->foff_est_hz = -fshift_max + foff_out; - fprintf(stderr,"syncing, shift %f,fine foff %f\n",fshift_max,-fshift_max + foff_out); - } - fprintf(stderr,"Best ratio %f freq %f offset %d fine %f\n",corr_ratio_max,fshift_max,frame_pos_est,ofdm->foff_est_hz); } } -// void ofdm_demod_coarse(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in){ -// const int SampsPerFrame = ofdm->config.SampsPerFrame; -// const int RxBufSize = ofdm->config.RxBufSize; -// const int Fs = ofdm->config.Fs; -// const int M = ofdm->config.M; -// const int Ncp = ofdm->config.Ncp; - -// int sync = ofdm->sync_count; -// int frame_point = ofdm->frame_point; -// int nin = ofdm->nin; -// int frame_pos_est; -// float fshift; -// float corr_ratio_max = 0; -// float corr_ratio; -// float foff_out; -// float fshift_max; -// int frame_pos_max = 0; - -// /* NCO for freq. shifting */ -// complex float shift_nco = 1; -// complex float shift_nco_dph = 0; - -// ofdm_demod(ofdm,rx_bits,rxbuf_in); -// frame_pos_est = coarse_sync(ofdm,ofdm->rxbuf,RxBufSize,&foff_out,&corr_ratio); -// frame_pos_est = (frame_pos_est % SampsPerFrame) - (M+Ncp) + 1; -// //fprintf(stderr,"frame_pos_est is %d\n",frame_pos_est); -// fprintf(stderr,"t_est is %d, foff_est is %d, ofdm_foff_est is %f\n",frame_pos_est,(int)foff_out,ofdm->foff_est_hz); - -// if(sync <= 0 ) { -// if(corr_ratio > 0.5){ -// if(frame_pos_est > 12){ -// ofdm->nin = frame_pos_est; -// } -// float delta_f = fabsf((float)foff_out - ofdm->foff_est_hz); -// if(delta_f > 3){ -// ofdm->foff_est_hz = (float)foff_out; -// } -// ofdm->sync_count = 3; -// } -// }else{ -// if(corr_ratio < 0.5){ -// ofdm->sync_count--; -// } -// float delta_f = fabsf((float)foff_out - ofdm->foff_est_hz); -// if(delta_f > 3){ -// ofdm->foff_est_hz = (float)foff_out; -// } -// } -// } - /* * ------------------------------------------ * ofdm_demod - Demodulates one frame of bits @@ -911,59 +432,36 @@ void ofdm_demod_coarse(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in){ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { complex float aphase_est_pilot_rect; - - const int Fs = ofdm->config.Fs; - const float Rs = ofdm->config.Rs; - const int Bps = ofdm->config.bps; - const int Nc = ofdm->config.Nc; - const int M = ofdm->config.M; - const int Ncp = ofdm->config.Ncp; - const int Ns = ofdm->config.Ns; - const int RowsPerFrame = ofdm->config.RowsPerFrame; - const int SampsPerFrame = ofdm->config.SampsPerFrame; - const int FtWindowWidth = ofdm->config.FtWindowWidth; - const int RxBufSize = ofdm->config.RxBufSize; - int nin = ofdm->nin; - - float aphase_est_pilot[Nc + 2]; - float aamp_est_pilot[Nc + 2]; + float aphase_est_pilot[OFDM_NC + 2]; + float aamp_est_pilot[OFDM_NC + 2]; float freq_err_hz; int i, j, k, rr, st, en, ft_est; - complex float * rxbuf = ofdm->rxbuf; - - - memmove(&rxbuf[0],&rxbuf[nin], (RxBufSize - nin) * sizeof(rxbuf[0])); - memcpy (&rxbuf[RxBufSize - nin],rxbuf_in, nin * sizeof(rxbuf[0])); - /* shift the buffer left based on nin */ - // for (i = 0, j = ofdm->nin; i < (RxBufSize - ofdm->nin); i++, j++) { - // ofdm->rxbuf[i] = ofdm->rxbuf[j]; - // } + for (i = 0, j = ofdm->nin; i < (OFDM_RXBUF - ofdm->nin); i++, j++) { + ofdm->rxbuf[i] = ofdm->rxbuf[j]; + } /* insert latest input samples onto tail of rxbuf */ - /* Note: COMP and complex float have the same memory layout */ - /* see iso/iec 9899:1999 sec 6.2.5.13 */ - // for (i = (RxBufSize - ofdm->nin), j = 0; i < RxBufSize; i++, j++) { - // ofdm->rxbuf[i] = rxbuf_in[j].real + rxbuf_in[j].imag * I; - // } + for (i = (OFDM_RXBUF - ofdm->nin), j = 0; i < OFDM_RXBUF; i++, j++) { + ofdm->rxbuf[i] = rxbuf_in[j].real + rxbuf_in[j].imag * I; + } /* * get user and calculated freq offset */ - float woff_est = TAU * ofdm->foff_est_hz / (float)(Fs); - float freq_offset = -1; - float ratio_out; + float woff_est = TAU * ofdm->foff_est_hz / OFDM_FS; + /* update timing estimate -------------------------------------------------- */ if (ofdm->timing_en == true) { /* update timing at start of every frame */ - st = ((M + Ncp) + SampsPerFrame) - floorf(FtWindowWidth / 2) + ofdm->timing_est; - en = st + SampsPerFrame - 1 + (M + Nc) + FtWindowWidth; + st = ((OFDM_M + OFDM_NCP) + OFDM_SAMPLESPERFRAME) - floorf(OFDM_FTWINDOWWIDTH / 2) + ofdm->timing_est; + en = st + OFDM_SAMPLESPERFRAME - 1 + (OFDM_M + OFDM_NCP) + OFDM_FTWINDOWWIDTH; complex float work[(en - st)]; @@ -973,41 +471,20 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { */ for (i = st, j = 0; i < en; i++, j++) { - /* TODO: Look into reducing cexpf usage */ work[j] = ofdm->rxbuf[i] * cexpf(-I * woff_est * i); } - // fprintf(stderr,"n>"); - ft_est = coarse_sync(ofdm, work, (en - st), &freq_offset,&ratio_out); - fprintf(stderr,"good ratio %f\n",ratio_out); - // fprintf(stderr,"\nx>"); - if(ratio_out < .7){ - ofdm->sync_count--; - } - //ofdm->foff_est_hz += (freq_offset/2); - ofdm->timing_est += (ft_est - ceilf(FtWindowWidth / 2)); + ft_est = coarse_sync(ofdm, work, (en - st)); + ofdm->timing_est += (ft_est - ceilf(OFDM_FTWINDOWWIDTH / 2)); if (ofdm->verbose > 1) { fprintf(stdout, " ft_est: %2d timing_est: %2d sample_point: %2d\n", ft_est, ofdm->timing_est, ofdm->sample_point); } - - //fprintf(stderr," freq_est: %f timing_est: %d sample_pt: %d\n",freq_offset,ofdm->timing_est,ofdm->sample_point); /* Black magic to keep sample_point inside cyclic prefix. Or something like that. */ - ofdm->sample_point = max(ofdm->timing_est + (Ncp / 4), ofdm->sample_point); - ofdm->sample_point = min(ofdm->timing_est + Ncp, ofdm->sample_point); - - - int ft_est_other = coarse_sync(ofdm,ofdm->rxbuf,RxBufSize,&freq_offset,&ratio_out); - ft_est_other = (ft_est_other % SampsPerFrame) - (M+Ncp) + 1; - // fprintf(stderr,"\n"); - //coarser_sync(ofdm, work, (en - st), &freq_offset, NULL); - - fprintf(stderr,"ratio is %f t_offset is %d f_off is %f freq_off is %f\n",ratio_out,ft_est_other,ofdm->foff_est_hz,freq_offset); - if(freq_offset > 1){ - //ofdm->foff_est_hz += freq_offset/4; - } + ofdm->sample_point = max(ofdm->timing_est + (OFDM_NCP / 4), ofdm->sample_point); + ofdm->sample_point = min(ofdm->timing_est + OFDM_NCP, ofdm->sample_point); } /* @@ -1049,8 +526,8 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { * The average of the four pilot symbols is our phase estimation. */ - for (i = 0; i < (Ns + 3); i++) { - for (j = 0; j < (Nc + 2); j++) { + for (i = 0; i < (OFDM_NS + 3); i++) { + for (j = 0; j < (OFDM_NC + 2); j++) { ofdm->rx_sym[i][j] = 0.0f + 0.0f * I; } } @@ -1059,10 +536,10 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { * "Previous" pilot symbol is one modem frame above. */ - st = (M + Ncp) + 1 + ofdm->sample_point; - en = st + M; + st = (OFDM_M + OFDM_NCP) + 1 + ofdm->sample_point; + en = st + OFDM_M; - complex float work[M]; + complex float work[OFDM_M]; /* down-convert at current timing instant---------------------------------- */ @@ -1098,9 +575,9 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { * In this routine we also process the current data symbols. */ - for (rr = 0; rr < (Ns + 1); rr++) { - st = (M + Ncp) + SampsPerFrame + (rr * (M + Ncp)) + 1 + ofdm->sample_point; - en = st + M; + for (rr = 0; rr < (OFDM_NS + 1); rr++) { + st = (OFDM_M + OFDM_NCP) + OFDM_SAMPLESPERFRAME + (rr * (OFDM_M + OFDM_NCP)) + 1 + ofdm->sample_point; + en = st + OFDM_M; /* down-convert at current timing instant---------------------------------- */ @@ -1147,8 +624,8 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { * We only want the "future" pilot symbol, to perform the averaging of all pilots. */ - st = (M + Ncp) + (3 * SampsPerFrame) + 1 + ofdm->sample_point; - en = st + M; + st = (OFDM_M + OFDM_NCP) + (3 * OFDM_SAMPLESPERFRAME) + 1 + ofdm->sample_point; + en = st + OFDM_M; /* down-convert at current timing instant---------------------------------- */ @@ -1167,7 +644,7 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { * +----------------------+ */ - dft(ofdm, ofdm->rx_sym[Ns + 2], work); + dft(ofdm, ofdm->rx_sym[OFDM_NS + 2], work); /* * We are finished now with the DFT and down conversion @@ -1185,19 +662,19 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { */ complex float freq_err_rect = conjf(vector_sum(ofdm->rx_sym[1], - Nc + 2)) * vector_sum(ofdm->rx_sym[Ns + 1], Nc + 2); + OFDM_NC + 2)) * vector_sum(ofdm->rx_sym[OFDM_NS + 1], OFDM_NC + 2); /* prevent instability in atan(im/re) when real part near 0 */ freq_err_rect = freq_err_rect + 1E-6f; - freq_err_hz = cargf(freq_err_rect) * Rs / (TAU * (float)Ns); + freq_err_hz = cargf(freq_err_rect) * OFDM_RS / (TAU * OFDM_NS); ofdm->foff_est_hz += (ofdm->foff_est_gain * freq_err_hz); } /* OK - now estimate and correct pilot phase ---------------------------------- */ - for (i = 0; i < (Nc + 2); i++) { + for (i = 0; i < (OFDM_NC + 2); i++) { aphase_est_pilot[i] = 10.0f; aamp_est_pilot[i] = 0.0f; } @@ -1208,7 +685,7 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { * Then average the phase surrounding each of the data symbols. */ - for (i = 1; i < (Nc + 1); i++) { + for (i = 1; i < (OFDM_NC + 1); i++) { complex float symbol[3]; for (j = (i - 1), k = 0; j < (i + 2); j++, k++) { @@ -1218,7 +695,7 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { aphase_est_pilot_rect = vector_sum(symbol, 3); for (j = (i - 1), k = 0; j < (i + 2); j++, k++) { - symbol[k] = ofdm->rx_sym[Ns + 1][j] * conjf(ofdm->pilots[j]); /* next pilot conjugate */ + symbol[k] = ofdm->rx_sym[OFDM_NS + 1][j] * conjf(ofdm->pilots[j]); /* next pilot conjugate */ } aphase_est_pilot_rect = aphase_est_pilot_rect + vector_sum(symbol, 3); @@ -1232,7 +709,7 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { aphase_est_pilot_rect = aphase_est_pilot_rect + vector_sum(symbol, 3); for (j = (i - 1), k = 0; j < (i + 2); j++, k++) { - symbol[k] = ofdm->rx_sym[Ns + 2][j] * ofdm->pilots[j]; /* last pilot */ + symbol[k] = ofdm->rx_sym[OFDM_NS + 2][j] * ofdm->pilots[j]; /* last pilot */ } aphase_est_pilot_rect = aphase_est_pilot_rect + vector_sum(symbol, 3); @@ -1250,18 +727,17 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { */ complex float rx_corr; - complex float rx_np[RowsPerFrame * Nc]; int abit[2]; int bit_index = 0; - for (rr = 0; rr < RowsPerFrame; rr++) { + for (rr = 0; rr < OFDM_ROWSPERFRAME; rr++) { /* * Note the i has an index of 1 to 16, so we ignore carriers 0 and 17. * * Also note we are using sym[2..8] or the seven data symbols. */ - for (i = 1; i < (Nc + 1); i++) { + for (i = 1; i < (OFDM_NC + 1); i++) { if (ofdm->phase_est_en == true) { rx_corr = ofdm->rx_sym[rr + 2][i] * cexpf(-I * aphase_est_pilot[i]); } else { @@ -1273,7 +749,7 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { * without the pilots. Thus, the name rx (no pilot) np */ - rx_np[(rr * Nc) + (i - 1)] = rx_corr; + ofdm->rx_np[(rr * OFDM_NC) + (i - 1)] = rx_corr; /* * Note even though amp ests are the same for each col, @@ -1281,18 +757,18 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { * so convenient to log them all */ - ofdm->rx_amp[(rr * Nc) + (i - 1)] = aamp_est_pilot[i]; + ofdm->rx_amp[(rr * OFDM_NC) + (i - 1)] = aamp_est_pilot[i]; /* * Note like amps in this implementation phase ests are the * same for each col, but we log them for each symbol anyway */ - ofdm->aphase_est_pilot_log[(rr * Nc) + (i - 1)] = aphase_est_pilot[i]; + ofdm->aphase_est_pilot_log[(rr * OFDM_NC) + (i - 1)] = aphase_est_pilot[i]; - if (Bps == 1) { + if (OFDM_BPS == 1) { rx_bits[bit_index++] = crealf(rx_corr) > 0.0f; - } else if (Bps == 2) { + } else if (OFDM_BPS == 2) { /* * Only one final task, decode what quadrant the phase * is in, and return the dibits @@ -1306,30 +782,21 @@ void ofdm_demod(struct OFDM *ofdm, int *rx_bits, COMP *rxbuf_in) { /* Adjust nin to take care of sample clock offset */ - ofdm->nin = SampsPerFrame; + ofdm->nin = OFDM_SAMPLESPERFRAME; if (ofdm->timing_en == true) { - int thresh = (M + Ncp) / 8; - int tshift = (M + Ncp) / 4; + int thresh = (OFDM_M + OFDM_NCP) / 8; + int tshift = (OFDM_M + OFDM_NCP) / 4; if (ofdm->timing_est > thresh) { - ofdm->nin = SampsPerFrame + tshift; + ofdm->nin = OFDM_SAMPLESPERFRAME + tshift; ofdm->timing_est -= tshift; ofdm->sample_point -= tshift; } else if (ofdm->timing_est < -thresh) { - ofdm->nin = SampsPerFrame - tshift; + ofdm->nin = OFDM_SAMPLESPERFRAME - tshift; ofdm->timing_est += tshift; ofdm->sample_point += tshift; } } - - int timing_est_log = ofdm->timing_est + 1; - int sample_point_log = ofdm->sample_point + 1; - - modem_probe_samp_f("foff_hz_log_c",&(ofdm->foff_est_hz),1); - modem_probe_samp_i("timing_est_log_c",&timing_est_log,1); - modem_probe_samp_i("sample_point_log_c",&sample_point_log,1); - modem_probe_samp_cft("rxbuf_log_c",ofdm->rxbuf,ofdm->config.RxBufSize); - modem_probe_samp_cft("rx_np_log_c",&rx_np[0],RowsPerFrame * Nc); } diff --git a/codec2-dev/src/ofdm_internal.h b/codec2-dev/src/ofdm_internal.h index bf1f1e83..1a1c576f 100644 --- a/codec2-dev/src/ofdm_internal.h +++ b/codec2-dev/src/ofdm_internal.h @@ -34,10 +34,8 @@ extern "C" { #include #include -#include #include "codec2_ofdm.h" -#include "kiss_fft.h" #ifndef M_PI #define M_PI 3.14159265358979323846f /* math constant */ @@ -45,55 +43,72 @@ extern "C" { #define TAU (2.0f * M_PI) /* mathematical constant */ +#define OFDM_NC 16 /* N Carriers */ +#define OFDM_TS 0.018f /* Symbol time */ +#define OFDM_RS (1.0f / OFDM_TS) /* Symbol rate */ +#define OFDM_FS 8000.0f /* Sample rate */ +#define OFDM_BPS 2 /* Bits per symbol */ +#define OFDM_TCP 0.002f /* ? */ +#define OFDM_NS 8 /* */ +#define OFDM_CENTRE 1500.0f /* Center frequency */ + +/* To prevent C99 warning */ + +#define OFDM_M 144 /* Samples per bare symbol (?) */ +#define OFDM_NCP 16 /* Samples per cyclic prefix */ + +#ifdef OLD_STYLE +/* This will produce a warning in C99 as (int) makes these variable */ + +#define OFDM_M ((int)(OFDM_FS / OFDM_RS)) +#define OFDM_NCP ((int)(OFDM_TCP * OFDM_FS)) +#endif + +/* ? */ +#define OFDM_FTWINDOWWIDTH 11 +/* Bits per frame (duh) */ +#define OFDM_BITSPERFRAME ((OFDM_NS - 1) * (OFDM_NC * OFDM_BPS)) +/* Rows per frame */ +#define OFDM_ROWSPERFRAME (OFDM_BITSPERFRAME / (OFDM_NC * OFDM_BPS)) +/* Samps per frame */ +#define OFDM_SAMPLESPERFRAME (OFDM_NS * (OFDM_M + OFDM_NCP)) + +#define OFDM_MAX_SAMPLESPERFRAME (OFDM_SAMPLESPERFRAME + (OFDM_M + OFDM_NCP)/4) +#define OFDM_RXBUF (3 * OFDM_SAMPLESPERFRAME + 3 * (OFDM_M + OFDM_NCP)) + + +/* Dummy struct for now, will contain constant configuration for OFDM modem */ struct OFDM_CONFIG{ - int32_t Nc; /* N carriers*/ - float Ts; /* Symbol time (S) (no CP) */ - float Rs; /* Symbol rate (S) (no CP) */ - int32_t Fs; /* Sample rate */ - int32_t bps; /* Bits per symbol */ - int32_t Tcp; /* Cyclic Prefix time (S) */ - int32_t Ns; /* Symbols per frame, including pilot */ - int32_t Fcenter; /* Center frequency */ - int32_t M; /* Samples per symbol (no CP) */ - int32_t Ncp; /* Samples of cyclic prefix */ - int32_t FtWindowWidth; /* ? */ - int32_t BitsPerFrame; /* Bits in a frame */ - int32_t SampsPerFrame; /* Samples in a frame */ - int32_t SampsPerFrameMax; /* Maximum samples per demod cycle */ - int32_t RxBufSize; /* RX buffer size */ - int32_t RowsPerFrame; /* Symbol depth per frame, no prefix */ + int a; }; -struct OFDM { - struct OFDM_CONFIG config; /* OFDM configuration (see above) */ - float foff_est_gain; /* ? */ - float foff_est_hz; /* Estimated frequency offset */ - - int verbose; /* Printing verbosity level */ - int sample_point; /* Fine timing offset */ - int timing_est; /* Ditto (?) (Coarse timing, maybe?) */ - int nin; /* Samples expected next demod cycle */ - int frame_point; - int sync_count; - - bool timing_en; /* Enable fine timing estimation */ - bool foff_est_en; /* Enable frequency offset estimation/tracking */ - bool phase_est_en; /* Enable Phase offset estimation */ - - complex float ** W; - complex float * pilots; - complex float * pilot_samples; - complex float * rxbuf; - float * w; - kiss_fft_cfg sync_fft_cfg; - - complex float * coarse_rxbuf; +struct OFDM { + struct OFDM_CONFIG config; + float foff_est_gain; + float foff_est_hz; + + int verbose; + int sample_point; + int timing_est; + int nin; + + bool timing_en; + bool foff_est_en; + bool phase_est_en; + + complex float pilot_samples[OFDM_M + OFDM_NCP]; + complex float W[OFDM_NC + 2][OFDM_M]; + complex float rxbuf[OFDM_RXBUF]; + complex float pilots[OFDM_NC + 2]; + float w[OFDM_NC + 2]; + /* Demodulator data */ - complex float ** rx_sym; - float * rx_amp; - float * aphase_est_pilot_log; + complex float rx_sym[OFDM_NS + 3][OFDM_NC + 2]; + complex float rx_np[OFDM_ROWSPERFRAME * OFDM_NC]; + float rx_amp[OFDM_ROWSPERFRAME * OFDM_NC]; + float aphase_est_pilot_log[OFDM_ROWSPERFRAME * OFDM_NC]; }; #ifdef __cplusplus diff --git a/codec2-dev/unittest/tofdm.c b/codec2-dev/unittest/tofdm.c index b886fef9..2ffff9a2 100644 --- a/codec2-dev/unittest/tofdm.c +++ b/codec2-dev/unittest/tofdm.c @@ -35,50 +35,15 @@ #include #include -#define MODEMPROBE_ENABLE - #include "ofdm_internal.h" #include "codec2_ofdm.h" #include "octave.h" #include "test_bits_ofdm.h" #include "comp_prim.h" -#include "modem_probe.h" - -#define OFDM_NC 16 - -#define NFRAMES 60 -//#define SAMPLE_CLOCK_OFFSET_PPM 100 -//#define FOFF_HZ 5.0f - - -#define OFDM_NCX 16 /* N Carriers */ -#define OFDM_TS 0.018f /* Symbol time */ -#define OFDM_RS (1.0f / OFDM_TS) /* Symbol rate */ -#define OFDM_FS 8000.0f /* Sample rate */ -#define OFDM_BPS 2 /* Bits per symbol */ -#define OFDM_TCP 0.002f /* ? */ -#define OFDM_NS 8 /* Symbols per frame (number of rows incl pilot) */ -#define OFDM_CENTRE 1500.0f /* Center frequency */ - -/* ? */ -#define OFDM_FTWINDOWWIDTH 11 -/* Bits per frame (duh) */ -#define OFDM_BITSPERFRAME ((OFDM_NS - 1) * (OFDM_NCX * OFDM_BPS)) -/* Rows per frame */ -#define OFDM_ROWSPERFRAME (OFDM_BITSPERFRAME / (OFDM_NCX * OFDM_BPS)) -/* Samps per frame */ -#define OFDM_SAMPLESPERFRAME (OFDM_NS * (OFDM_M + OFDM_NCP)) - -#define OFDM_MAX_SAMPLESPERFRAME (OFDM_SAMPLESPERFRAME + (OFDM_M + OFDM_NCP)/4) -#define OFDM_RXBUF (3 * OFDM_SAMPLESPERFRAME + 3 * (OFDM_M + OFDM_NCP)) - - - -/* To prevent C99 warning */ - -#define OFDM_M 144 /* Samples per bare symbol (?) */ -#define OFDM_NCP 16 /* Samples per cyclic prefix */ +#define NFRAMES 30 +#define SAMPLE_CLOCK_OFFSET_PPM 100 +#define FOFF_HZ 0.5f /*---------------------------------------------------------------------------*\ @@ -152,12 +117,10 @@ static void freq_shift(COMP rx_fdm_fcorr[], COMP rx_fdm[], float foff, COMP *fof int main(int argc, char *argv[]) { + int samples_per_frame = ofdm_get_samples_per_frame(); + int max_samples_per_frame = ofdm_get_max_samples_per_frame(); - - struct OFDM *ofdm = ofdm_create(OFDM_CONFIG_700D); - int samples_per_frame = ofdm_get_samples_per_frame(ofdm); - int max_samples_per_frame = ofdm_get_max_samples_per_frame(ofdm); - + struct OFDM *ofdm; COMP tx[samples_per_frame]; /* one frame of tx samples */ int rx_bits[OFDM_BITSPERFRAME]; /* one frame of rx bits */ @@ -166,8 +129,7 @@ int main(int argc, char *argv[]) int tx_bits_log[OFDM_BITSPERFRAME*NFRAMES]; COMP tx_log[samples_per_frame*NFRAMES]; - //COMP rx_log[samples_per_frame*NFRAMES]; - COMP * rx_log; + COMP rx_log[samples_per_frame*NFRAMES]; COMP rxbuf_in_log[max_samples_per_frame*NFRAMES]; COMP rxbuf_log[OFDM_RXBUF*NFRAMES]; COMP rx_sym_log[(OFDM_NS + 3)*NFRAMES][OFDM_NC + 2]; @@ -182,10 +144,7 @@ int main(int argc, char *argv[]) FILE *fout; int f,i,j; - int sample_clock_offset_ppm = 100; - float foff_hz = 0.1; - int nframes = 30; - + ofdm = ofdm_create(OFDM_CONFIG_700D); assert(ofdm != NULL); /* Main Loop ---------------------------------------------------------------------*/ @@ -203,30 +162,19 @@ int main(int argc, char *argv[]) /* tx vector logging */ - memcpy(&tx_bits_log[OFDM_BITSPERFRAME*f], test_bits_ofdm, sizeof(int)*OFDM_BITSPERFRAME); - memcpy(&tx_log[samples_per_frame*f], tx, sizeof(COMP)*samples_per_frame); + memcpy(&tx_bits_log[OFDM_BITSPERFRAME*f], test_bits_ofdm, sizeof(int)*OFDM_BITSPERFRAME); + memcpy(&tx_log[samples_per_frame*f], tx, sizeof(COMP)*samples_per_frame); } /* --------------------------------------------------------*\ Channel \*---------------------------------------------------------*/ - /*fs_offset(rx_log, tx_log, samples_per_frame*NFRAMES, sample_clock_offset_ppm); + fs_offset(rx_log, tx_log, samples_per_frame*NFRAMES, SAMPLE_CLOCK_OFFSET_PPM); COMP foff_phase_rect = {1.0f, 0.0f}; - freq_shift(rx_log, rx_log, foff_hz, &foff_phase_rect, samples_per_frame * NFRAMES);*/ - - FILE *cplin = fopen("tofdm_rx_vec","r"); - fseek(cplin,0,SEEK_END); - size_t file_len = ftell(cplin); - rewind(cplin); - size_t rx_samps_in = file_len/sizeof(COMP); - rx_log = malloc(rx_samps_in*sizeof(COMP)); - fread(rx_log,sizeof(COMP),rx_samps_in,cplin); - fclose(cplin); - //fprintf(stderr,"Read %d\n",s); - + freq_shift(rx_log, rx_log, FOFF_HZ, &foff_phase_rect, samples_per_frame * NFRAMES); /* --------------------------------------------------------*\ Demod @@ -235,108 +183,96 @@ int main(int argc, char *argv[]) /* Init rx with ideal timing so we can test with timing estimation disabled */ int Nsam = samples_per_frame*NFRAMES; - Nsam = rx_samps_in; int prx = 0; - //int nin = samples_per_frame + 2*(OFDM_M+OFDM_NCP); - int nin = ofdm_get_nin(ofdm); + int nin = samples_per_frame + 2*(OFDM_M+OFDM_NCP); + int lnew; COMP rxbuf_in[max_samples_per_frame]; - /*for (i=0; irxbuf[OFDM_RXBUF-nin+i] = rx_log[prx].real + I*rx_log[prx].imag; - }*/ + } int nin_tot = 0; - modem_probe_init("ofdm","ofdm_test.txt"); - /* disable estimators for initial testing */ ofdm_set_verbose(ofdm, true); ofdm_set_timing_enable(ofdm, true); ofdm_set_foff_est_enable(ofdm, true); ofdm_set_phase_est_enable(ofdm, true); - //ofdm->foff_est_hz = 10; - - bool cont = true; - int ptr = 0; - for(f=0; cont; f++) { - // /* For initial testng, timing est is off, so nin is always - // fixed. TODO: we need a constant for rxbuf_in[] size that - // is the maximum possible nin */ - - // nin = ofdm_get_nin(ofdm); - - // nin_tot += nin; - - // assert(nin <= max_samples_per_frame); - - // /* Insert samples at end of buffer, set to zero if no samples - // available to disable phase estimation on future pilots on - // last frame of simulation. */ - - // if ((Nsam-prx) < nin) { - // lnew = Nsam-prx; - // } else { - // lnew = nin; - // } - // //printf("nin: %d prx: %d lnew: %d\n", nin, prx, lnew); - // for(i=0; i Nsam){ - // cont = false; - // break; - // } - //assert(prx <= max_samples_per_frame*NFRAMES); - - //fprintf(stderr,"NIN:%d\n",nin); + + for(f=0; f= Nsam){ - break; - } + assert(nin <= max_samples_per_frame); - memset(rxbuf_in,0,sizeof(COMP) * nin); + /* Insert samples at end of buffer, set to zero if no samples + available to disable phase estimation on future pilots on + last frame of simulation. */ - for(i=0; i= Nsam){ - cont = false; - //break; + if (lnew) { + for(i=0; irx_sym[i][j]); - rx_sym_log[(OFDM_NS + 3)*f+i][j].imag = cimagf(ofdm->rx_sym[i][j]); - } + for(i=0; irxbuf[i]); + rxbuf_log[OFDM_RXBUF*f+i].imag = cimagf(ofdm->rxbuf[i]); + } + + for (i = 0; i < (OFDM_NS + 3); i++) { + for (j = 0; j < (OFDM_NC + 2); j++) { + rx_sym_log[(OFDM_NS + 3)*f+i][j].real = crealf(ofdm->rx_sym[i][j]); + rx_sym_log[(OFDM_NS + 3)*f+i][j].imag = cimagf(ofdm->rx_sym[i][j]); } } + + /* note corrected phase (rx no phase) is one big linear array for frame */ + + for (i = 0; i < OFDM_ROWSPERFRAME*OFDM_NC; i++) { + rx_np_log[OFDM_ROWSPERFRAME*OFDM_NC*f + i].real = crealf(ofdm->rx_np[i]); + rx_np_log[OFDM_ROWSPERFRAME*OFDM_NC*f + i].imag = cimagf(ofdm->rx_np[i]); + } + /* note phase/amp ests the same for each col, but check them all anyway */ - if(f < NFRAMES){ - for (i = 0; i < OFDM_ROWSPERFRAME; i++) { - for (j = 0; j < OFDM_NC; j++) { - phase_est_pilot_log[OFDM_ROWSPERFRAME*f+i][j] = ofdm->aphase_est_pilot_log[OFDM_NC*i+j]; - rx_amp_log[OFDM_ROWSPERFRAME*OFDM_NC*f+OFDM_NC*i+j] = ofdm->rx_amp[OFDM_NC*i+j]; - } + + for (i = 0; i < OFDM_ROWSPERFRAME; i++) { + for (j = 0; j < OFDM_NC; j++) { + phase_est_pilot_log[OFDM_ROWSPERFRAME*f+i][j] = ofdm->aphase_est_pilot_log[OFDM_NC*i+j]; + rx_amp_log[OFDM_ROWSPERFRAME*OFDM_NC*f+OFDM_NC*i+j] = ofdm->rx_amp[OFDM_NC*i+j]; } } - modem_probe_samp_i("rx_bits_log_c",rx_bits,OFDM_BITSPERFRAME); + + foff_hz_log[f] = ofdm->foff_est_hz; + timing_est_log[f] = ofdm->timing_est + 1; /* offset by 1 to match Octave */ + sample_point_log[f] = ofdm->sample_point + 1; /* offset by 1 to match Octave */ + + memcpy(&rx_bits_log[OFDM_BITSPERFRAME*f], rx_bits, sizeof(rx_bits)); } /*---------------------------------------------------------*\ @@ -344,9 +280,6 @@ int main(int argc, char *argv[]) by tofdm.m Octave script \*---------------------------------------------------------*/ - modem_probe_close(); - - fout = fopen("tofdm_out.txt","wt"); assert(fout != NULL); fprintf(fout, "# Created by tofdm.c\n"); @@ -354,10 +287,16 @@ int main(int argc, char *argv[]) octave_save_int(fout, "tx_bits_log_c", tx_bits_log, 1, OFDM_BITSPERFRAME*NFRAMES); octave_save_complex(fout, "tx_log_c", (COMP*)tx_log, 1, samples_per_frame*NFRAMES, samples_per_frame*NFRAMES); octave_save_complex(fout, "rx_log_c", (COMP*)rx_log, 1, samples_per_frame*NFRAMES, samples_per_frame*NFRAMES); + octave_save_complex(fout, "rxbuf_in_log_c", (COMP*)rxbuf_in_log, 1, nin_tot, nin_tot); + octave_save_complex(fout, "rxbuf_log_c", (COMP*)rxbuf_log, 1, OFDM_RXBUF*NFRAMES, OFDM_RXBUF*NFRAMES); octave_save_complex(fout, "rx_sym_log_c", (COMP*)rx_sym_log, (OFDM_NS + 3)*NFRAMES, OFDM_NC + 2, OFDM_NC + 2); octave_save_float(fout, "phase_est_pilot_log_c", (float*)phase_est_pilot_log, OFDM_ROWSPERFRAME*NFRAMES, OFDM_NC, OFDM_NC); octave_save_float(fout, "rx_amp_log_c", (float*)rx_amp_log, 1, OFDM_ROWSPERFRAME*OFDM_NC*NFRAMES, OFDM_ROWSPERFRAME*OFDM_NC*NFRAMES); - + octave_save_float(fout, "foff_hz_log_c", foff_hz_log, NFRAMES, 1, 1); + octave_save_int(fout, "timing_est_log_c", timing_est_log, NFRAMES, 1); + octave_save_int(fout, "sample_point_log_c", sample_point_log, NFRAMES, 1); + octave_save_complex(fout, "rx_np_log_c", (COMP*)rx_np_log, 1, OFDM_ROWSPERFRAME*OFDM_NC*NFRAMES, OFDM_ROWSPERFRAME*OFDM_NC*NFRAMES); + octave_save_int(fout, "rx_bits_log_c", rx_bits_log, 1, OFDM_BITSPERFRAME*NFRAMES); fclose(fout); ofdm_destroy(ofdm); -- 2.25.1