fft base fast xcorr is done. see fft_corr.m

Author: sdr-x

Matlab/CellSearch.m

@@ -57,7 +57,7 @@
     disp(['Input averaged abs: ' num2str( mean(abs([real(capbuf_pbch) imag(capbuf_pbch)])) )]);
 
     disp('sampling_ppm_f_search_set_by_pss: try ... ... ');
-    [period_ppm, dynamic_f_search_set, xc, ~, ~, ~, ~, extra_info] = sampling_ppm_f_search_set_by_pss(capbuf_pbch.', f_search_set, pss_fo_set, sampling_carrier_twist, pss_peak_max_reserve, num_pss_period_try, combined_pss_peak_range, par_th, num_peak_th);
+    [period_ppm, dynamic_f_search_set, xc, ~, ~, ~, ~, extra_info] = sampling_ppm_f_search_set_by_pss(capbuf_pbch.', f_search_set, td_pss, pss_fo_set, sampling_carrier_twist, pss_peak_max_reserve, num_pss_period_try, combined_pss_peak_range, par_th, num_peak_th);
     if sampling_carrier_twist==0
         if period_ppm == inf
             disp('No valid PSS is found at pre-proc phase! Please try again.');

Matlab/fft_corr.m

@@ -0,0 +1,34 @@
+function [corr_store, fo_search_set_new] = fft_corr(s, td_pss, fo_search_set)
+% % s column vectors -- 153600
+% % td_pss column vectors
+% % fo_search_set row vectors
+
+fo_search_set_new = fo_search_set; % it happens to be equal, because 1.92e6/153600 = 12.5Hz under 5kHz step size
+sampling_rate = 1.92e6;
+
+len = length(s);
+
+freq_step = sampling_rate/len;
+
+len_pss = size(td_pss, 1);
+num_pss = size(td_pss, 2);
+num_fo = length(fo_search_set);
+num_fo_pss = num_fo*num_pss;
+% len_short = len - (len_pss-1);
+
+corr_store = zeros(len, num_fo_pss);
+
+s = fft(s);
+fd_pss = fft(conj(td_pss(end:-1:1,:)./len_pss), len, 1);
+
+for i = 1 : num_fo_pss
+    pss_idx = floor((i-1)/num_fo) + 1;
+    fo_idx = i - (pss_idx-1)*num_fo;
+    fo = fo_search_set(fo_idx);
+    fd_fo_shift_len = fo/freq_step;
+    corr_store(:, i) = s.*circshift(fd_pss(:,pss_idx), [fd_fo_shift_len,0]);
+%     corr_store(:, i) = abs( ifft(s.*circshift(fd_pss(:,pss_idx), [fd_fo_shift_len,0]), [], 1) ).^2;
+%     corr_store(:, i) = real( corr_store(:, i).*conj(corr_store(:, i)) );
+end
+corr_store = abs( ifft(corr_store, len, 1) ).^2;
+corr_store = corr_store(len_pss:end,:);

Matlab/sampling_ppm_f_search_set_by_pss.m

@@ -2,13 +2,13 @@
 % Find out LTE PSS in the signal stream and correct sampling&carrier error.
 % A script of project: https://github.com/JiaoXianjun/rtl-sdr-LTE
 
-function [ppm, f_set, xc, fo_idx_set, pss_idx_set, fo_pss_idx_set, fo_with_all_pss_idx, extra_info] = sampling_ppm_f_search_set_by_pss(s, fo_search_set, pss_fo_set, sampling_carrier_twist, max_reserve, num_pss_period_try, combined_pss_peak_range, par_th, num_peak_th)
+function [ppm, f_set, xc, fo_idx_set, pss_idx_set, fo_pss_idx_set, fo_with_all_pss_idx, extra_info] = sampling_ppm_f_search_set_by_pss(s, fo_search_set, td_pss, pss_fo_set, sampling_carrier_twist, max_reserve, num_pss_period_try, combined_pss_peak_range, par_th, num_peak_th)
 % sampling period PPM! not sampling frequency PPM!
 
 % fo_search_set = -100e3 : 5e3 : 100e3; % -100kHz ~ 100 kHz with 5kHz step size
 % pss_fo_set = pss_fo_set_gen(td_pss, fo_search_set);
 
-combined_pss_peak_range_half = floor(combined_pss_peak_range/2);
+% combined_pss_peak_range_half = floor(combined_pss_peak_range/2);
 extra_info = [];
 
 len_pss = size(pss_fo_set, 1);
@@ -17,13 +17,21 @@
 len = length(s);
 len_short = len - (len_pss-1);
 
-corr_store = zeros(len_short, num_fo_pss);
+% corr_store = zeros(len_short, num_fo_pss);
+% 
+% tic;
+% for i=1:num_fo_pss
+%     tmp_corr = abs( filter(pss_fo_set(end:-1:1, i), 1, s) ).^2;
+%     tmp_corr = tmp_corr(len_pss:end);
+%     corr_store(:,i) = tmp_corr;
+% end
+% cost_time1 = toc;
+% disp(['PSS xcorr cost ' num2str(cost_time1)]);
 
-for i=1:num_fo_pss
-    tmp_corr = abs( filter(pss_fo_set(end:-1:1, i), 1, s) ).^2;
-    tmp_corr = tmp_corr(len_pss:end);
-    corr_store(:,i) = tmp_corr;
-end
+tic;
+[corr_store, fo_search_set] = fft_corr(s, td_pss, fo_search_set);
+cost_time2 = toc;
+disp(['PSS xcorr cost ' num2str(cost_time2)]);
 
 if sampling_carrier_twist==1
     ppm = inf;

Matlab/test_CellSearch.m

@@ -31,9 +31,9 @@
 
 cell_info = cell(1, length(test_source_info));
 for i = test_sp : test_ep
-    if isempty( strfind(test_source_info(i).filename, 'dimitri') )
-        continue;
-    end
+%     if isempty( strfind(test_source_info(i).filename, 'dimitri') )
+%         continue;
+%     end
     disp(test_source_info(i).filename);
     coef_pbch = pbch_filter_coef_gen(test_source_info(i).fs);
 

Matlab/test_fft_corr.m

@@ -0,0 +1,43 @@
+% function test_fft_corr
+clear all;
+close all;
+[~, td_pss] = pss_gen;
+td_pss = td_pss(:,1).';
+len_pss = length(td_pss);
+
+fo_search_set = 5e3;
+
+sampling_rate = 1.92e6;
+td_pss_fo = td_pss.*exp(1i.*2.*pi.*(1./sampling_rate).*(0:(len_pss-1)).*fo_search_set);
+
+conj_td_pss_fo = conj(td_pss_fo);
+
+len = 153600;
+s = randn(1, len) + 1i.*randn(1, len);
+
+freq_step = sampling_rate/len;
+
+shift_len = fo_search_set/freq_step;
+
+fft_conj_td_pss_fo = fft(conj_td_pss_fo, len);
+
+fft_conj_td_pss = fft(conj(td_pss), len);
+
+plot(angle(fft_conj_td_pss_fo)); hold on;
+plot(angle(circshift(fft_conj_td_pss, [0, -shift_len])), 'r');
+
+fft_s = fft(s);
+
+% a = filter(conj_td_pss_fo(end:-1:1), 1, s);
+% a = conv(s, conj_td_pss_fo);
+% a(1: len_pss-1) = a(1: len_pss-1) + a(len+1 : end);
+a = filter(conj_td_pss_fo, 1, s);
+
+b = ifft( fft_s.*circshift(fft_conj_td_pss, [0, -shift_len]) );
+% b = ifft( fft_s.*fft_conj_td_pss_fo );
+
+% b = b(len_pss:end);
+
+figure;
+plot(angle(a)); hold on;
+plot(angle(b),'r');