diff --git a/ofdm/ofdm_rx.m b/ofdm/ofdm_rx.m index a08bc28..6a5b65e 100644 --- a/ofdm/ofdm_rx.m +++ b/ofdm/ofdm_rx.m @@ -1,9 +1,9 @@ -function ofdm_rx(filename, foff, flip_spec, SNR, with_cir, show_H) +function ofdm_rx(filename, foff, flip_spec, SNR, show_H) % -% [sym_out, H, W] = ofdm_rx(filename, foff, flip_spec, SNR, with_cir, show_H) -% Example: ofdm_rx('drm-15435-if.wav', -174.65, 0, 400, 0, 1) -% Example: ofdm_rx('xv.wav', 0, 0, 400, 0, 1) -% or ofdm_rx('yv.wav', 0, 0, 400, 0, 1) +% [sym_out, H, W] = ofdm_rx(filename, foff, flip_spec, SNR, show_H) +% Example: ofdm_rx('drm-15435-if.wav', -174.65, 0, 400, 1) +% Example: ofdm_rx('xv.wav', 0, 0, 400, 1) +% or ofdm_rx('yv.wav', 0, 0, 400, 1) close all; @@ -12,204 +12,68 @@ drm_bw = '10k'; [ofdm_drm_params, ofdm_spec_occ] = drm_params(drm_mode, drm_bw); -pre_mix = 0; fs = 12000; N = 256; Ng = 64; + [carrier_pilot, phi_pilot, mag_pilot] = getRefFreq(drm_mode); [carrier_time, phi_time, mag_time] = getRefTime(drm_mode); -h_cir = [0 0 .1 -.2 .5 .2 -.1 0 0 ]; +% Read wav file +[fsize fs4] = wavread(filename, 'size') +[xwav, fs4, nbits] = wavread(filename, 0); +wav_numCh = fsize(2); +wav_size = fsize(1); -% Generate signal -[xwav fs_ nbits] = wavread(filename); - -[len ndim] = size(xwav); - -if (ndim == 2) - fprintf('Mode: Baseband\n'); - x = xwav(:, 1) + j*xwav(:, 2); - x = x.*exp(-j*(2*pi*(foff/fs_)*(0:len-1)')); -elseif (ndim == 1) +isModePassband = 0; +if (wav_numCh == 1) fprintf('Mode: Passband\n'); - x = xwav - mean(xwav); - x = xwav.*exp(-j*(2*pi*(0.25-foff/fs_)*(0:len-1)')); - hlp = FIRCalcLowpass(12000/fs_, 15); - x = filter(hlp, 1, x); - - % Downsample - x = x(1:4:length(x)); + omega_xwav = (0.25-foff/fs4); + isModePassband = 1; +elseif (wav_numCh == 2) + fprintf('Mode: Baseband\n'); + omega_xwav = (foff/fs4); else - error('Cannot read file!'); + error('Cannot read file'); end -if with_cir - x = filter(h_cir, 1, x); -end +buf_recv = buffer(); +buf_Z = buffer(); -x = x - mean(x); -if flip_spec - x = conj(x); -end -x = x./max(abs(x)); -xout = [real(x) imag(x)]; -wavwrite(xout,fs,nbits,'ofdm_rx.wav'); +fileChunkSize = N; +fileCunkCounter = 0; +phi_xwav = 0; +filter_h = FIRCalcLowpass(0.25, 15); +filter_z = []; -% Modulate on carrier -xup = upsample(x, 4); -hlp = FIRCalcLowpass(0.25, 129); -xup = filter(hlp, 1, xup); +fileRemain = wav_size; -y = exp(j*2*pi.*(0:length(xup)-1)'*0.25) .* xup; +% Timing coarse stuff +timingCoarse = struct('Z', 0, 'sync_counter', 20, 'sync', 0, 'cp_off', 0, 'cp_off_', [], 'Nd', 0); -wavwrite(y,4*fs,nbits,'ofdm_rx_passband.wav'); +% CFO coarse integer stuff +cfoCoarseInteger = struct('Z', 0, 'sync_counter', 20, 'sync', 0, 'df', 0); -x = awgn(x, SNR, 'measured'); -K = fix((length(x))/(N+Ng))-1; +% CFO coarse fractional stuff +cfoCoarseFractional = struct('Z', 0, 'sync_counter', 200, 'sync', 0, 'err', 0, 'err_', [], 'df', 0, 'df_', []); -% Timing coarse -fprintf('Coarse timing acquisition\n'); -cp_off = 0; -cp_off_ = []; -Z = 0; -for k=1:20 - xp = x((k-1)*(N+Ng)+1:k*(N+Ng)); - [cp_off, Nd, valid, Z] = timing_coarse(xp, N, Ng, Ng, Z); - if (valid) - cp_off_ = [cp_off_ cp_off]; - end -end; -cp_off = round(median(cp_off_)) -Nd = Nd +% CFO Fine stuff +cfoFine = struct('Z', 0, 'sync_counter', 200, 'sync', 0, 'err', 0, 'err_', [], 'df', 0, 'df_', [], 'Xlast', [], 'phi', 0, 'bin_track', 1); -% Integer CFO coarse -fprintf('Coarse integer CFO acquisition\n'); -fdet_coarse_Z = 0; -Z_cfa_pre = 0; -df_acq_coarse_int_ = []; -valid_latency_count = 20; -for k=1:200 - xp = x((k-1)*(N+Ng)+cp_off+1:k*(N+Ng)+cp_off); - [d_bin, valid, fdet_coarse_Z] = cfo_coarse_int(xp(Ng+1:N+Ng), N, carrier_pilot, 11, fdet_coarse_Z); - if valid - df_acq_coarse_int = -d_bin(1)/N*fs; - if valid_latency_count > 0 - valid_latency_count = valid_latency_count - 1; - else - break; - end - end -end +% Frame timing stuff +frameTiming = struct('Z', 0, 'sync_counter', 20, 'sync', 0, 'kc_max', -1000, 'n_max', 0); -fprintf('Coarse fractional CFO acquisition\n'); -% CFO coarse -df_acq_coarse_fract_ = []; -for k=1:200 - xp = x((k-1)*(N+Ng)+cp_off+1:k*(N+Ng)+cp_off); - [domega_f, Z_cfa_pre] = cfo_coarse_fract(xp, N, Ng, Ng, Nd, Z_cfa_pre); - df_acq_coarse_fract_ = [df_acq_coarse_fract_ domega_f*fs]; -end; -df_acq_c = df_acq_coarse_int + mean(df_acq_coarse_fract_(10:200)) -df_acq_c_ = df_acq_coarse_int + df_acq_coarse_fract_; - -if pre_mix - x1 = x.*exp(j*(2*pi*df_acq_c/fs.*(0:length(x)-1)')); - df_acq_fine = 0; -else - x1 = x; - df_acq_fine = df_acq_c; -end - -% CFO Fine -fprintf('Fine fractional CFO acquisition\n'); -Xlast = []; -df_acq_fine_err = 0; -df_acq_fine_err_ = []; -df_acq_fine_ = []; -kn = fs/(2*pi*(1+Ng/N))/N; -x_ = []; -Hp_ = []; -Z_ = []; -phi = 0.0; -bin_track = 1; -sym_pilot = exp(j*2*pi*phi_pilot'/1024); - -bin_pilot = c2i(N, carrier_pilot); -for k=1:K, - xp = x1((k-1)*(N+Ng)+cp_off+1:k*(N+Ng)+cp_off); - xp = xp.*exp(j*(2*pi*df_acq_fine/fs*(0:N+Ng-1)' + phi)); - phi = phi + 2*pi*df_acq_fine/fs*(N+Ng); - xp = xp(Ng+1:N+Ng); - x_ = [x_' xp']'; - X = fft(xp); - if (~isempty(Xlast)) - Z = Xlast .* conj(X); - df_acq_fine_err = kn*angle(Z(bin_pilot)); - df_acq_fine_err_ = [df_acq_fine_err_ df_acq_fine_err]; - phi_p = angle(X(bin_pilot)); - mag_p = abs(X(bin_pilot)); - df_acq_fine = df_acq_fine + 0.1*mean(df_acq_fine_err(bin_track)); - df_acq_fine_ = [df_acq_fine_ df_acq_fine]; - Hp = X(bin_pilot)./conj(sym_pilot); - Hp_ = [Hp_ Hp]; - Z_ = [Z_ X]; - end - Xlast = X; -end - -figure; -plot(0:length(cp_off_)-1, cp_off_, 0:length(cp_off_)-1, cp_off*ones(size(cp_off_)), 'r-'); grid; title('cp_{off}') - -df_acq_fine_ = df_acq_fine_'; - -figure; -subplot(3, 1, 1) -plot(0:length(df_acq_c_)-1, df_acq_c_); grid; title('df_{coarse} (Acquisition)') - -subplot(3, 1, 2) -plot(0:length(df_acq_fine_)-1, df_acq_fine_); grid; title('df_{fine} (Acquisition)') - -subplot(3, 1, 3) -plot(0:length(df_acq_fine_err_)-1, df_acq_fine_err_); grid; title('Error(df_{fine}) (Acquisition)') - -figure; -Z_pilot = Z_(bin_pilot, 200:length(Z_))'; -plot(real(Z_pilot), imag(Z_pilot), '+'); grid; title('SymPilot_{Freq}'); -maxZ = max(max(abs(Z_pilot))); -if (maxZ < 1) - axis([-1 1 -1 1]); -else - axis([-maxZ maxZ -maxZ maxZ]); -end - -fprintf('Frame timing acquisition\n'); -%% Detect start of DRM frame -bins = c2i(N, carrier_time); -kc_max = -1000; -n_max = 0; -skip_syms = ofdm_drm_params.nspf*4; -for n=0:ofdm_drm_params.nspf-1 - frame_timing_Z = 0; - k = skip_syms+n; - for m=1:2 - [kc valid frame_timing_Z] = timing_frame(Z_(bins, k)', bins, ofdm_drm_params.nspf, frame_timing_Z); - k = k + ofdm_drm_params.nspf; - end - if valid - if kc > kc_max - n_max = n; - kc_max = kc; - end - end -end -frame_start = n_max - -fprintf('Channel estimation and symbol reception\n'); +% Equalization stuff +MIN_ABS_H = 1e-10; +equalizer = struct('Z', 0, 'sync_counter', 20, 'sync', 0, 's', 0, 'm', 0, 'kk', 1, 'eps_abs_h', MIN_ABS_H); [W_syms, W_pilots] = calcWiener(drm_mode, drm_bw); -figure; -% Detect scattered gain pilots +H_syms = []; +H_pilots = []; +H_ = []; +df_track_ = []; + symbols_per_frame = ofdm_drm_params.nspf; symbols_to_delay = ofdm_drm_params.y; frames_per_window = 2*ofdm_drm_params.y; @@ -235,100 +99,263 @@ fac_cell_list = ... [23, 35, 47, 65, 77], [] }; -s = 0; -m = 0; -H_syms = []; -H_pilots = []; -kk = 1; -H_ = []; -df_track_ = []; -MIN_ABS_H = 1e-10; -eps_abs_h = MIN_ABS_H; +while(fileRemain > fileChunkSize) + + n0 = fileCunkCounter*fileChunkSize+1; + n1 = (fileCunkCounter+1)*fileChunkSize; + [xwav] = wavread(filename, [n0 n1]); + + if (~isModePassband) + xwav = xwav(:, 1) + j*xwav(:, 2); + end + xwav = xwav - mean(xwav); -start = skip_syms+frame_start; + x = xwav.*exp(-j*(2*pi*omega_xwav*(0:fileChunkSize-1)' + phi_xwav)); + phi_xwav = mod(phi_xwav + 2*pi*omega_xwav*fileChunkSize, 2*pi); + + if flip_spec + x = conj(x); + end -for k=start:length(Z_) - - %shifted symbol index - n = rem(m+symbols_to_delay-1, ofdm_drm_params.y); - - [carrier_gain, phi_gain, mag_gain] = getRefGain(ofdm_drm_params, ofdm_spec_occ, s); - carrier_gain_bins = c2i(N, carrier_gain); - H_ = transpose(Z_(carrier_gain_bins, k))./(mag_gain.*exp(j*2*pi*phi_gain/1024)); - - H_raw(carrier_gain_bins) = (H_); - H_gain_cell{m+1} = (H_); - - str = sprintf('Frame %d, s=%d, m=%d, n=%d', k, s, m, n); - if kk > frames_per_window - ii = 1+mod(m+symbols_to_delay-1+(0:frames_per_window-1), frames_per_window); - - H_gain = [H_gain_cell{ii}]; - H_syms = H_gain*W_syms{n+1}; - H_pilots = (H_gain)*W_pilots{n+1}; - - % Get the correctly delayed symbols - sym = Z_(carrier_indexes, k-symbols_to_delay).'; - - % Get the correctly delayed raw channel estimate - Hr = H_gain_cell{mod(m-1, frames_per_window)+1}; - - % get frequency error - % correlate wiener filtered channel with raw estimate - df_track = fs/(2*pi*N)*angle(H_pilots * Hr' + eps_abs_h); - df_track_ = [df_track_ df_track]; - - % check for too small values - rms_H = sqrt( abs( H_syms*H_syms' )/length(H_syms) ); - eps_abs_h = max( rms_H/20, MIN_ABS_H ); %quick and dirty fixed settings - H_too_small_index = find( abs(H_syms) 0) + timingCoarse.sync_counter = timingCoarse.sync_counter - 1; + else + timingCoarse.sync = 1; end - hold off; - - pause(0.01); - else - fac_c = fac_cell_list{mod(s-symbols_to_delay, symbols_per_frame)+1}; - if (~isempty(fac_c)) - fac_i = fac_c - ofdm_spec_occ.kmin + 1; - hold on; - plot(sym_eq(fac_i), '+', 'MarkerSize', 4); - pause(0.01); - grid; title('FAC'); axis([-1 1 -1 1]); - hold off; + end; + if (timingCoarse.sync) + timingCoarse.cp_off = round(median(timingCoarse.cp_off_)); + timingCoarse.Nd = Nd; + Nd = timingCoarse.Nd + cp_off = timingCoarse.cp_off + + % Reset buffer to Cp offset + buf_recv.setReadIndex(timingCoarse.cp_off); + end + end + + % Integer CFO coarse + if (timingCoarse.sync) && (~cfoCoarseInteger.sync) + x = buf_recv.read(N+Ng); + if (~isempty(x)) + [d_bin, valid, cfoCoarseInteger.Z] = cfo_coarse_int(x(Ng+1:N+Ng), N, carrier_pilot, 11, cfoCoarseInteger.Z); + if valid + cfoCoarseInteger.df = -d_bin(1)/N*fs; + if (cfoCoarseInteger.sync_counter > 0) + cfoCoarseInteger.sync_counter = cfoCoarseInteger.sync_counter - 1; + else + cfoCoarseInteger.sync = 1; + end end end + end + + % Fractional CFO coarse + if (cfoCoarseInteger.sync) && (~cfoCoarseFractional.sync) + x = buf_recv.read(N+Ng); + if (~isempty(x)) + [err, cfoCoarseFractional.Z] = cfo_coarse_fract(x, N, Ng, Ng, Nd, cfoCoarseFractional.Z); + cfoCoarseFractional.err_ = [cfoCoarseFractional.err_ err*fs]; + + if (cfoCoarseFractional.sync_counter > 0) + cfoCoarseFractional.sync_counter = cfoCoarseFractional.sync_counter - 1; + else + cfoCoarseFractional.sync = 1; + cfoCoarseFractional.df = cfoCoarseInteger.df + mean(cfoCoarseFractional.err_); + cfoCoarseFractional.df_ = cfoCoarseInteger.df + cfoCoarseFractional.err_; + cfoFine.df = cfoCoarseFractional.df; + end end - end - kk = kk + 1; - m = mod(m+1, frames_per_window); - s = mod(s + 1, symbols_per_frame); + end + + % CFO Fine + if (cfoCoarseFractional.sync) + kn = fs/(2*pi*(1+Ng/N))/N; + sym_pilot = exp(j*2*pi*phi_pilot'/1024); + bin_pilot = c2i(N, carrier_pilot); + + x = buf_recv.read(N+Ng); + if (~isempty(x)) + x = x.*exp(j*(2*pi*cfoFine.df/fs*(0:N+Ng-1)' + cfoFine.phi)); + cfoFine.phi = cfoFine.phi + 2*pi*cfoFine.df/fs*(N+Ng); + X = fft(x(Ng+1:N+Ng)); + if (~isempty(cfoFine.Xlast)) + Z = cfoFine.Xlast .* conj(X); + cfoFine.err = kn*angle(Z(bin_pilot)); + cfoFine.err_ = [cfoFine.err_ cfoFine.err]; + cfoFine.df = cfoFine.df + 0.1*mean(cfoFine.err(cfoFine.bin_track)); + cfoFine.df_ = [cfoFine.df_ cfoFine.df]; + buf_Z.write(X); + end + cfoFine.Xlast = X; + end + + bins = c2i(N, carrier_time); + skip_syms = ofdm_drm_params.nspf*4; + % Frame timing acquisition + if (~frameTiming.sync) && (buf_Z.len() >= (N*(ofdm_drm_params.nspf+skip_syms))) + buf_Z.setReadIndex(N*skip_syms); + for n=0:ofdm_drm_params.nspf-1 + frameTiming.Z = 0; + k = skip_syms+n; + Z = buf_Z.read(N); + if ~isempty(Z) + for m=1:2 + [kc valid frameTiming.Z] = timing_frame(Z(bins)', bins, ofdm_drm_params.nspf, frameTiming.Z); + k = k + ofdm_drm_params.nspf; + end + if valid + if kc > frameTiming.kc_max + frameTiming.n_max = n; + frameTiming.kc_max = kc; + end + end + end + end + if ~isempty(Z) + buf_Z.setReadIndex(N*(frameTiming.n_max)); + frameTiming.sync = 1; + frame_start = frameTiming.n_max + end + end + + if (frameTiming.sync) + + Z = transpose(buf_Z.read(N)); + if ~isempty(Z) + + equalizer.n = rem(equalizer.m + symbols_to_delay-1, ofdm_drm_params.y); + + [carrier_gain, phi_gain, mag_gain] = getRefGain(ofdm_drm_params, ofdm_spec_occ, equalizer.s); + carrier_gain_bins = c2i(N, carrier_gain); + H_ = Z(carrier_gain_bins)./(mag_gain.*exp(j*2*pi*phi_gain/1024)); + + H_raw(carrier_gain_bins) = (H_); + H_gain_cell{equalizer.m+1} = (H_); + + str = sprintf('Frame %d, s=%d, m=%d, n=%d', equalizer.kk, equalizer.s, equalizer.m, equalizer.n); + + if equalizer.kk > frames_per_window + + ii = 1+mod(equalizer.m+symbols_to_delay-1+(0:frames_per_window-1), frames_per_window); + + H_gain = [H_gain_cell{ii}]; + H_syms = H_gain*W_syms{equalizer.n+1}; + H_pilots = (H_gain)*W_pilots{equalizer.n+1}; + + % Get the correctly delayed symbols + Z = buf_Z.readAt(N, -N*symbols_to_delay); + sym = Z(carrier_indexes).'; + + % Get the correctly delayed raw channel estimate + Hr = H_gain_cell{mod(equalizer.m-1, frames_per_window)+1}; + + % get frequency error + % correlate wiener filtered channel with raw estimate + df_track = fs/(2*pi*N)*angle(H_pilots * Hr' + equalizer.eps_abs_h); + df_track_ = [df_track_ df_track]; + + % check for too small values + rms_H = sqrt( abs( H_syms*H_syms' )/length(H_syms) ); + equalizer.eps_abs_h = max( rms_H/20, MIN_ABS_H ); %quick and dirty fixed settings + H_too_small_index = find( abs(H_syms)