From 90231d49ae31ae8158bd6651a0efcd7ab05fd9cc Mon Sep 17 00:00:00 2001 From: yeokwanggoo-stack Date: Tue, 29 Jul 2025 21:45:44 +0900 Subject: [PATCH] Initial Commit. --- Functions/00. General/cal_doppler_shift.m | 22 + Functions/00. General/cal_isodoppler.m | 46 ++ Functions/00. General/cal_winloss.m | 46 ++ Functions/00. General/findpeak2D.m | 50 ++ Functions/00. General/fn_add_data_tip.m | 8 + Functions/00. General/gen_freqgrid.m | 49 ++ Functions/00. General/peaks2.m | 239 +++++++ Functions/00. General/poly2_fit.m | 8 + .../cal_mutual_coherence.m | 24 + .../01. Direction finding/cal_sv_ambi_func.m | 28 + Functions/01. Direction finding/df_1D_MUSIC.m | 44 ++ Functions/01. Direction finding/df_FFT.m | 40 ++ .../01. Direction finding/df_FFT2d_mar510.m | 62 ++ .../01. Direction finding/df_FFT_mar510.m | 60 ++ Functions/01. Direction finding/df_OMP.m | 61 ++ .../df_Phase_monopulse.m | 5 + Functions/01. Direction finding/df_cbf.m | 67 ++ .../01. Direction finding/gen_sig_for_df.m | 124 ++++ .../01. Direction finding/gen_virarray.m | 36 + Functions/03. FMCW radar design/SNR_cal.m | 68 ++ .../03. FMCW radar design/azi_ant_pat.mat | Bin 0 -> 1003 bytes Functions/03. FMCW radar design/beat2rng.m | 11 + Functions/03. FMCW radar design/bw2rngsep.m | 11 + Functions/03. FMCW radar design/cal_nf.m | 38 + .../03. FMCW radar design/cal_rdcoupling.m | 11 + .../03. FMCW radar design/cal_rtmcoupling.m | 7 + Functions/03. FMCW radar design/dechirping.m | 10 + Functions/03. FMCW radar design/dop2spd.m | 7 + .../03. FMCW radar design/elev_ant_pat.mat | Bin 0 -> 1064 bytes .../03. FMCW radar design/gen_fastramp_fmcw.m | 68 ++ .../gen_fastramp_fmcw2.m | 71 ++ Functions/03. FMCW radar design/rng2beat.m | 15 + Functions/03. FMCW radar design/rng2time.m | 11 + Functions/03. FMCW radar design/rngres2bw.m | 11 + Functions/03. FMCW radar design/spd2dop.m | 7 + .../cal_rts_attenuation.m | 10 + Functions/04. Target modelling/gen_tgt.m | 10 + System_Analysis_v2_0.m | 666 ++++++++++++++++++ 38 files changed, 2051 insertions(+) create mode 100644 Functions/00. General/cal_doppler_shift.m create mode 100644 Functions/00. General/cal_isodoppler.m create mode 100644 Functions/00. General/cal_winloss.m create mode 100644 Functions/00. General/findpeak2D.m create mode 100644 Functions/00. General/fn_add_data_tip.m create mode 100644 Functions/00. General/gen_freqgrid.m create mode 100644 Functions/00. General/peaks2.m create mode 100644 Functions/00. General/poly2_fit.m create mode 100644 Functions/01. Direction finding/cal_mutual_coherence.m create mode 100644 Functions/01. Direction finding/cal_sv_ambi_func.m create mode 100644 Functions/01. Direction finding/df_1D_MUSIC.m create mode 100644 Functions/01. Direction finding/df_FFT.m create mode 100644 Functions/01. Direction finding/df_FFT2d_mar510.m create mode 100644 Functions/01. Direction finding/df_FFT_mar510.m create mode 100644 Functions/01. Direction finding/df_OMP.m create mode 100644 Functions/01. Direction finding/df_Phase_monopulse.m create mode 100644 Functions/01. Direction finding/df_cbf.m create mode 100644 Functions/01. Direction finding/gen_sig_for_df.m create mode 100644 Functions/01. Direction finding/gen_virarray.m create mode 100644 Functions/03. FMCW radar design/SNR_cal.m create mode 100644 Functions/03. FMCW radar design/azi_ant_pat.mat create mode 100644 Functions/03. FMCW radar design/beat2rng.m create mode 100644 Functions/03. FMCW radar design/bw2rngsep.m create mode 100644 Functions/03. FMCW radar design/cal_nf.m create mode 100644 Functions/03. FMCW radar design/cal_rdcoupling.m create mode 100644 Functions/03. FMCW radar design/cal_rtmcoupling.m create mode 100644 Functions/03. FMCW radar design/dechirping.m create mode 100644 Functions/03. FMCW radar design/dop2spd.m create mode 100644 Functions/03. FMCW radar design/elev_ant_pat.mat create mode 100644 Functions/03. FMCW radar design/gen_fastramp_fmcw.m create mode 100644 Functions/03. FMCW radar design/gen_fastramp_fmcw2.m create mode 100644 Functions/03. FMCW radar design/rng2beat.m create mode 100644 Functions/03. FMCW radar design/rng2time.m create mode 100644 Functions/03. FMCW radar design/rngres2bw.m create mode 100644 Functions/03. FMCW radar design/spd2dop.m create mode 100644 Functions/04. Target modelling/cal_rts_attenuation.m create mode 100644 Functions/04. Target modelling/gen_tgt.m create mode 100644 System_Analysis_v2_0.m diff --git a/Functions/00. General/cal_doppler_shift.m b/Functions/00. General/cal_doppler_shift.m new file mode 100644 index 0000000..09fe5cb --- /dev/null +++ b/Functions/00. General/cal_doppler_shift.m @@ -0,0 +1,22 @@ +function fd = cal_doppler_shift(va, vb, ra, rb, fc) + % va, vb : 각각 물체 a, b의 속도 벡터 [vx vy vz] + % ra, rb : 각각 물체 a, b의 위치 벡터 [x y z] + % fc : 송신 주파수 (Hz) + + c = 3e8; % 빛의 속도 (m/s) + + % 단위 방향 벡터 r_hat 계산 + r_ab = rb - ra; % a에서 b로 향하는 벡터 + r_hat = r_ab / norm(r_ab); % 단위 벡터 + + % 속도 차이 벡터 계산 + v_rel = vb - va; + + % 도플러 시프트 계산 + fd = (fc / c) * dot(v_rel, r_hat); +end + + + + + diff --git a/Functions/00. General/cal_isodoppler.m b/Functions/00. General/cal_isodoppler.m new file mode 100644 index 0000000..65e04ff --- /dev/null +++ b/Functions/00. General/cal_isodoppler.m @@ -0,0 +1,46 @@ +clear all +close all +clc + +azi = -90 : 0.1 : 90; +elv = -15 : 0.1 : 15; + +% Reference point (origin) +x_rdr = 0; +y_rdr = 0; +z_rdr = 0; + +vx_rdr = 0; +vy_rdr = 70/3.6; +vz_rdr = 0; + +x_tgt_str = 4; +y_tgt_str = -10; +z_tgt_str = 0; + +x_tgt_mov = 0; +y_tgt_mov = -10; +z_tgt_mov = 0; + +vx_tgt_mov = 0; +vy_tgt_mov = 1.3905; +vz_tgt_mov = 0; + + +pos_rdr = [x_rdr, y_rdr, z_rdr]; +vel_rdr = [vx_rdr, vy_rdr, vz_rdr]; + +pos_tgt_str = [x_tgt_str, y_tgt_str, z_tgt_str]; +pos_tgt_mov = [x_tgt_mov, y_tgt_mov, z_tgt_mov]; +vel_tgt_mov = [vx_tgt_mov, vy_tgt_mov, vz_tgt_mov]; +vel_tgt_str = [0 0 0]; + +rel_r_vec_str = pos_tgt_str - pos_rdr +rel_r_vec_mov = pos_tgt_mov - pos_rdr +rel_v_vec_mov = vel_tgt_mov - vel_rdr +rel_v_vec_str = vel_tgt_str - vel_rdr + +rng_str = norm(rel_r_vec_str) +rng_mov = norm(rel_r_vec_mov) +rel_vel_mov = rel_v_vec_mov * rel_r_vec_mov'/norm(rel_r_vec_mov) +rel_vel_str = rel_v_vec_str * rel_r_vec_str'/norm(rel_r_vec_str) \ No newline at end of file diff --git a/Functions/00. General/cal_winloss.m b/Functions/00. General/cal_winloss.m new file mode 100644 index 0000000..c128f52 --- /dev/null +++ b/Functions/00. General/cal_winloss.m @@ -0,0 +1,46 @@ +function [worst_scalloping_loss_dB, avg_scalloping_loss_dB, SNR_loss_dB] = cal_winloss(win) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculating losses of windowing +% Start : 23.08.25 +% End : 23.08.25 +% developed by Kwanggoo Yeo +% +% Description +% - Input +% - 1) win [vector] : window function +% +% - Output +% - 1) worst_scalloping_loss_dB [scalar], [dB] : scalloping loss in worst case +% - 2) avg_scalloping_loss_dB [scalar], [dB] : scallopoing loss in average case [uniform distribution] +% - 3) SNR_loss_dB [scalar], [dB] : SNR processing gain loss +% +% +% History +% (23.08.25) Completed +% +% Referece +% - 1) Mark A. Richards, "Fundamentals of Radar Signal Processing 1st edition", p.257 +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +if size(win, 1) > 1 + win = win.'; +end + +N = length(win); + +worst_scalloping_loss_dB = mag2db( abs( sum( win .* exp(-1i * pi / N * [0:N-1]))) / sum(win)); + +SNR_loss_dB = pow2db( sum(win)^2 / (N * sum(win.^2)) ); + +findex = linspace(-pi/N, pi/N, 10000); + +for idx = 1 : length(findex) + val(idx) = (abs( sum( win .* exp(-1i * findex(idx) * [0:N-1]))) / sum(win)); +end +avg_scalloping_loss_dB = mag2db(sum(val)/10000); + +end + diff --git a/Functions/00. General/findpeak2D.m b/Functions/00. General/findpeak2D.m new file mode 100644 index 0000000..b830d0a --- /dev/null +++ b/Functions/00. General/findpeak2D.m @@ -0,0 +1,50 @@ +function [peakindex2D, peakval2D] = findpeak2D(data2D, th) + +[Nrow, Ncol] = size(data2D); + +peakindex2D = zeros(Nrow, Ncol); +peakval2D = zeros(Nrow, Ncol); + +if isempty(th) + for rowidx = 1 : Nrow + for colidx = 1 : Ncol + indexset = [rowidx-1 colidx; rowidx+1 colidx; rowidx colidx-1; rowidx colidx+1]; + det_indexset = (indexset(:,1) > 0) .* (indexset(:,1) < Nrow+1) .* (indexset(:,2) > 0) .* (indexset(:,2) < Ncol+1); + det_indexset = find(det_indexset > 0); + cond_peak = 0; + index_iter_len = length(det_indexset); + for indexset_idx = 1 : length(det_indexset) + if data2D(rowidx, colidx) >= data2D(indexset(det_indexset(indexset_idx),1), indexset(det_indexset(indexset_idx),2)) + cond_peak = cond_peak + 1; + end + end + + if cond_peak == index_iter_len + peakindex2D(rowidx, colidx) = 1; + peakval2D(rowidx, colidx) = data2D(rowidx, colidx); + end + end + end +else + for rowidx = 1 : Nrow + for colidx = 1 : Ncol + indexset = [rowidx-1 colidx; rowidx+1 colidx; rowidx colidx-1; rowidx colidx+1]; + det_indexset = (indexset(:,1) > 0) .* (indexset(:,1) < Nrow+1) .* (indexset(:,2) > 0) .* (indexset(:,2) < Ncol+1); + det_indexset = find(det_indexset > 0); + cond_peak = 0; + index_iter_len = length(det_indexset); + for indexset_idx = 1 : length(det_indexset) + if data2D(rowidx, colidx) >= data2D(indexset(det_indexset(indexset_idx),1), indexset(det_indexset(indexset_idx),2)) + cond_peak = cond_peak + 1; + end + end + + if cond_peak == index_iter_len && data2D(rowidx, colidx) >= th + peakindex2D(rowidx, colidx) = 1; + peakval2D(rowidx, colidx) = data2D(rowidx, colidx); + end + end + end +end + +end \ No newline at end of file diff --git a/Functions/00. General/fn_add_data_tip.m b/Functions/00. General/fn_add_data_tip.m new file mode 100644 index 0000000..fb50f78 --- /dev/null +++ b/Functions/00. General/fn_add_data_tip.m @@ -0,0 +1,8 @@ +function fn_add_data_tip(plot_data, name_of_data_tip, data, index) + row = dataTipTextRow(name_of_data_tip, data); + if ~isa(row.Value, 'double') + row.Value = string(row.Value); + end + plot_data.DataTipTemplate.DataTipRows(index) = row; +end + diff --git a/Functions/00. General/gen_freqgrid.m b/Functions/00. General/gen_freqgrid.m new file mode 100644 index 0000000..9dfdae2 --- /dev/null +++ b/Functions/00. General/gen_freqgrid.m @@ -0,0 +1,49 @@ +function freq_grid = gen_freqgrid(N, Fs, opt) +% gen_freqgrid Generate frequency grid +% freq_grid = gen_freqgrid(N, Fs, opt) generate an N point +% frequency grid according to sample rate Fs. This grid matches +% the operation used in fftshift(opt = 1) or not(opt = 0). +% +% % Example: +% % Create a 16 point frequency grid for a sample rate of 10 +% % Hz used in fftshift. +% +% freq_grid = gen_freqgrid(16, 10, 1) + +% set 'CenterDC' in psdfreqvec to true preserves Nyquist point, +% which does not match our processing to the data since we use +% fftshift. + +% adopted from psdfreqvec + +%% Checking 'opt' in input + if ~(opt == 0 || opt == 1) + disp('Option must be 0(No fftshift) or 1(fftshift)'); + return; + end + +%% Generating freqeuncy grid + % freq_grid = fftshift(psdfreqvec(... + % 'Npts',N,'Fs',Fs,'Range','whole')); + freq_res = Fs/N; + freq_grid = (0:N-1).'*freq_res; + + if opt == 1 + % linspace(freq_offset-fs/2, freq_offset+fs/2*(fft_len-2)/fft_len, fft_len); + Nyq = Fs/2; + half_res = freq_res/2; + if rem(N,2) % odd + idx = 1:(N-1)/2; + halfpts = (N+1)/2; + freq_grid(halfpts) = Nyq-half_res; + freq_grid(halfpts+1) = Nyq+half_res; + else + idx = 1:N/2; + hafpts = N/2+1; + freq_grid(hafpts) = Nyq; + end + freq_grid(N) = Fs-freq_res; + freq_grid = fftshift(freq_grid); + freq_grid(idx) = freq_grid(idx)-Fs; + end +end \ No newline at end of file diff --git a/Functions/00. General/peaks2.m b/Functions/00. General/peaks2.m new file mode 100644 index 0000000..c3b6b33 --- /dev/null +++ b/Functions/00. General/peaks2.m @@ -0,0 +1,239 @@ +function [pks,locs_y,locs_x]=peaks2(data,varargin) +% Find local peaks in 2D data. +% Syntax chosen to be as close as possible to the original Matlab +% 'findpeaks' function but not require any additional toolbox. +% +% SYNTAX: +% pks=peaks(data) finds local peaks. +% +% [pks,locs_y,locs_x]=peaks(data) finds local peaks and their array +% coordinates. +% +% [pks,locs_y,locs_x]=peaks(...,'MinPeakHeight',{scalar value}) only retains +% those peaks which are equal to or greater than this absolute value. +% +% [pks,locs_y,locs_x]=peaks(...,'Threshold',{scalar value}) only retains +% those peaks that are higher than their immediate surroundings by this value. +% +% [pks,locs_y,locs_x]=peaks(...,'MinPeakDistance',{scalar value}) finds peaks +% separated by more than the specified minimum CARTESIAN peak distance (a +% circle around the peak). It starts from the strongest peak and goes +% iteratively lower. Any peak 'shadowed' in the vicinity of a stronger +% peak is discarded. +% +% ALGORITHM: +% A peak is considered to be a data point strictly greater than its +% immediate neighbors. You can change this condition to 'greater or equal' +% in the code, but be aware that in such case, it might create false +% detections in flat areas, but these can be accounted for by introduction +% of a small Threshold value. +% +% Even though this function is shared here free for use and any +% modifications you might find useful, I would appreciate if you would +% quote me in case you are going to use this function for any non-personal +% tasks. +% (C) Kristupas Tikuisis 2023. + + +%% Initial data check +% Let's simplify the function. Let it work on 1D or 2D data only, and check +% if the input data satisfies this criteria: +if ~ismatrix(data) + error('Only 1D (vectors) or 2D (matrices) data accepted.'); +end + + +%% Locate all peaks +% A peak is a data point HIGHER than its immediate neighbors. There are 8 +% around eaxh point, and we will go through each. Oh yes, no escaping that. +% +% To introduce as little of intermediate variables and keep their memory +% footprint as low as possible, I will introduce 2 logical arrays: one to +% mark the peaks and be iteratively updated until we check all its +% neighbors; and another temporal variable just to prepare data for +% comparison (mind about the edge points which do not have any neighbors!): +ispeak=false(size(data)); % to store peak flags. +isgreater=true(size(data)); % to store comparison result for one particular neighbor. + +% Now start analyzing every data point. +% +% 1st neighbor immediatelly to the left: +ispeak=([true(size(data,1),1) [data(:,2:end)>data(:,1:end-1)]]); % for this case, we can update the peak array directly. +% +% 2nd neighbor at the top-left: +isgreater(2:end,2:end)=(data(2:end,2:end)>data(1:end-1,1:end-1)); % this time, due to points on the diagonal, a temporary array will have to be involved. +ispeak=ispeak&isgreater; % time to update the peak array. +% +% 3rd neighbor immediatelly at the top: +ispeak=ispeak&([true(1,size(data,2)); (data(2:end,:)>data(1:end-1,:))]); % once again, for this case, we can update the peak array directly. +% +% 4th neighbor at the top-right: +isgreater=true(size(data)); % rebuild a fresh array for comparison... +isgreater(2:end,1:end-1)=(data(2:end,1:end-1)>data(1:end-1,2:end)); +ispeak=ispeak&isgreater; +% +% 5th neighbor immediatelly to the right: +ispeak=ispeak&([(data(:,1:end-1)>data(:,2:end)) true(size(data,1),1)]); % once again, for this case, we can update the peak array directly. +% +% 6th neighbor to the bottom-right: +isgreater=true(size(data)); +isgreater(1:end-1,1:end-1)=(data(1:end-1,1:end-1)>data(2:end,2:end)); +ispeak=ispeak&isgreater; +% +% 7th neighbor immediatelly at the bottom: +ispeak=ispeak&([(data(1:end-1,:)>data(2:end,:)); true(1,size(data,2))]); % once again, for this case, we can update the peak array directly. +% +% 8th neighbor at the bottom-left: +isgreater=true(size(data)); +isgreater(1:end-1,2:end)=(data(1:end-1,2:end)>data(2:end,1:end-1)); +ispeak=ispeak&isgreater; +% +% Discard the temporary variable: +clear isgreater +% By now, raw peak indentification is completed. + + +%% Return final results +% First, essential raw peak location steps. Peak values and locations: +locs=find(ispeak); %clear ispeak +pks=data(locs); +% Note that LINEAR indices have been returned. Let's turn them to array +% indices for final output: +[locs_y,locs_x]=ind2sub(size(data),locs); + + +%% Perform customary post-processing determined by optional function parameters. +% First, check if any parameters were supplied: +if isempty(varargin) + return +end +% ...then a quality check - the parameters should come in pairs, therefore +% the length should be even: +if mod(length(varargin),2)~=0 + warning('Optional name-value parameters should come in pairs. Something is missing. Peak search will proceed with default values.'); + return +end +% ...after this step, we can sort the input into names (parameters) and +% values: +params=varargin(1:2:end); +vals=varargin(2:2:end); +% ...last quality check - all even values should be CHAR entries specifying +% a parameter to be adjusted: +ischarparam=cellfun(@ischar,params); +if any(~ischarparam) + warning('Some parameters are not written as characters and not recognisable. They should always come is name(char)-value pairs. Peak search will proceed with default values.'); + return +end; clear ischarparam varargin + +%-------------------------------------------------------------------------- +% No go over all supplied parameters and check the found peaks accordingly. + +% 1. MinPeakHeight - absolute minimum value for a peak: +isparm=find(cellfun(@(x)isequal(x,'MinPeakHeight'),params),1); +if ~isempty(isparm) + % Locate which peaks satisfy this condition: + suitable=(pks>=vals{isparm}); + % ...and only keep those: + pks=pks(suitable); + locs=locs(suitable); + clear suitable +end + +% 2. Threshold - peak must be greater than its neighbours by this value. +isparm=find(cellfun(@(x)isequal(x,'Threshold'),params),1); +if ~isempty(isparm) + % For this, we will need to convert linear indices to array indices: + [row_y,col_x]=ind2sub(size(data),locs); + % These will be the original indices. + + % This is how array coordinates would change relatively around each + % peak (y,x): + % (-1,-1) (-1,0) (-1,+1) + % ( 0,-1) ( 0,0) ( 0,+1) + % (+1,-1) (+1,0) (+1,+1) + % ...turned into vectors disregarding the (0,0), the center data point: + delta_y=[-1 -1 -1 0 0 +1 +1 +1]; + delta_x=[-1 0 +1 -1 +1 -1 0 +1]; + % Let's add these deltas to the detected peak positions to get the + % coordinates of their immediate neighbors: + neighbor_locs_y=row_y+delta_y; + neighbor_locs_x=col_x+delta_x; clear row_y col_x delta_x delta_y + % ...don't forget to check for unrealistic indices beyond array + % borders: + neighbor_locs_y(neighbor_locs_y<1)=1; + neighbor_locs_y(neighbor_locs_y>size(data,1))=size(data,1); + neighbor_locs_x(neighbor_locs_x<1)=1; + neighbor_locs_x(neighbor_locs_x>size(data,2))=size(data,2); + % ...convert to linear indices: + neighbor_locs=sub2ind(size(data),neighbor_locs_y,neighbor_locs_x); + clear neighbor_locs_y neighbor_locs_x + + % So we have neighbor values by now. Are our peaks higher than those by + % the set Threshold value? + suitable=(data(locs)-vals{isparm}>=data(neighbor_locs)); + % Now check for those cases when by mistake (earlier step for checking + % for indices beyond array boundaries) a peak itself is taken as a + % neighbor as well: + suitable(data(neighbor_locs)==data(locs))=true; + + % Only keep those is they are greater than ALL neighbors (in other + % words, those elements where NONE are lesser): + suitable=~any(~suitable,2); clear neighbor_locs + + % Final step - locate suitable element indices: + suitable=find(suitable); + + % That's it, modify the output array: + pks=pks(suitable); + locs=locs(suitable); clear suitable +end + +% 3. 'MinPeakDistance' +isparm=find(cellfun(@(x)isequal(x,'MinPeakDistance'),params),1); +if ~isempty(isparm) + + % First, sort the peaks in order of amplitude: + [pks_sorted,idx]=sort(pks,'descend'); + locs_sorted=locs(idx); clear idx + + % The flow is as follows: start from the highest peak and discard any + % other peaks closer than the CARTESIAN MinPeakDistance (that is, the + % CIRCLE around the peak is going to be checked); then continue until + % the whole list (updated iteratively as items might get removed) has + % been checked. + + % Convert locations to array indices: + [row_y,col_x]=ind2sub(size(data),locs_sorted); + + % Start from the highest peak: + this_peak=1; + while this_peak<(length(pks_sorted)+1) + + % Cartesian distances to ALL its remaining & yet unchecked neighbors + % (including itself): + dist=sqrt((row_y-row_y(this_peak)).^2+(col_x-col_x(this_peak)).^2); + % Now simply check which neighbors are WITHIN the MinPeakDistance + % BUT also nonzero (the peak should not be compared to itself): + within=( (dist<=vals{isparm}) & (dist~=0) ); + % ...and delete those entries satisfying the condition: + pks_sorted(within)=[]; + locs_sorted(within)=[]; + row_y(within)=[]; col_x(within)=[]; + + % Update the peak counter: + this_peak=this_peak+1; + + end; clear this_peak within dist row_y col_x + + % Update the peak location and value list: + pks=pks_sorted; + locs=locs_sorted; +end + +% I haven't figured out how to do it more elegantly, but turn the default +% linear indices to array indices: +[locs_y,locs_x]=ind2sub(size(data),locs); + +%========================================================================== +% End of the entire function +end \ No newline at end of file diff --git a/Functions/00. General/poly2_fit.m b/Functions/00. General/poly2_fit.m new file mode 100644 index 0000000..2f7f0ab --- /dev/null +++ b/Functions/00. General/poly2_fit.m @@ -0,0 +1,8 @@ +function [del_x, itp_y] = poly2_fit(y_1, y0, y1) + +del_x = -1/2 * ((y1) - (y_1)) / ((y_1) - 2 * (y0) + (y_1)); + +itp_y = 1/2 * ((del_x-1) * del_x * (y_1) - 2 * (del_x - 1) * (del_x + 1) * (y0) + (del_x+1) * del_x * (y_1)); + +end + diff --git a/Functions/01. Direction finding/cal_mutual_coherence.m b/Functions/01. Direction finding/cal_mutual_coherence.m new file mode 100644 index 0000000..9193479 --- /dev/null +++ b/Functions/01. Direction finding/cal_mutual_coherence.m @@ -0,0 +1,24 @@ +function [mc_val] = cal_mutual_coherence(dic_mat, plot_flag) + +[~, N] = size(dic_mat); + +for idx = 1 : N + for jdx = 1 : N + if idx ~= jdx + mc(idx, jdx) = abs(dic_mat(:, idx)' * dic_mat(:, jdx)) / norm(dic_mat(:, idx))^2; + else + mc(idx, jdx) = 0; + end + end +end + +if plot_flag == 1 + figure + imagesc(mc) +end + +mc_val = max(max(mc)); + + +end + diff --git a/Functions/01. Direction finding/cal_sv_ambi_func.m b/Functions/01. Direction finding/cal_sv_ambi_func.m new file mode 100644 index 0000000..6ab1f77 --- /dev/null +++ b/Functions/01. Direction finding/cal_sv_ambi_func.m @@ -0,0 +1,28 @@ +function [sv_ambi_val] = cal_sv_ambi_func(sv_1, sv_2) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculting ambiguity(normalized correlation) of two input steering vectors +% Start : 23.08.25 +% End : 23.08.25 +% developed by Kwanggoo Yeo +% +% Description +% - Input +% - 1) sv_1 [vector], [-] : the first steering vector fc +% - 2) sv_2 [vector], [-] : the second steering vector fc +% +% - Output +% - 1) sv_ambi_val [scalar], [-] : The ambiguity(normalized correlation) of two input steering vectors +% +% History +% (23.08.25) Completed +% +% Referece +% - +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +sv_ambi_val = abs(sv_1' * sv_2) / (norm(sv_1) * norm(sv_2)); + +end + diff --git a/Functions/01. Direction finding/df_1D_MUSIC.m b/Functions/01. Direction finding/df_1D_MUSIC.m new file mode 100644 index 0000000..1deb292 --- /dev/null +++ b/Functions/01. Direction finding/df_1D_MUSIC.m @@ -0,0 +1,44 @@ +function [esti_ang_deg, P_MUSIC, test_ang_deg] = df_1D_MUSIC(snapshot, Ntarget, array_struct, fov_deg, unit_ang_deg , plot_flag) + +% Number of Snapshots +Nsnap = size(snapshot, 2); + +% Estimated Correlation matrix of X +Rx = (1/Nsnap) * (snapshot * snapshot'); + +% MUSIC algorithm + % Eigen decomposition + [Q, ~, ~] = svd(Rx); + + % Noise subspace + Q_n = Q(:, Ntarget+1:end); + + % MUSIC algorithm + test_ang_deg = -fov_deg : unit_ang_deg : fov_deg - unit_ang_deg; + P_MUSIC = zeros(1, length(test_ang_deg)); + + if size(array_struct.eff_ch_loc, 1) == 1 + array_struct.eff_ch_loc = array_struct.eff_ch_loc.'; + end + + Q_n_square = Q_n * Q_n'; + for ang_idx = 1 : length(test_ang_deg) + z = exp(1i * 2 * pi / array_struct.lambda_c * array_struct.unit * -sind(test_ang_deg(ang_idx))); + test_sv = z.^array_struct.eff_ch_loc; + P_MUSIC(ang_idx) = 1 / (test_sv' * Q_n_square * test_sv); + end + +if plot_flag == 1 + figure() + plot(test_ang_deg, pow2db(abs(P_MUSIC))); + grid on + xlabel('Test angle [deg]'); + ylabel('Magnitude [dB]'); +end + +[pks, pks_ang, ~, ~] = findpeaks(pow2db(abs(P_MUSIC)), test_ang_deg); +[~, order] = sort(pks, 'descend'); +esti_ang_deg = pks_ang(order(1:Ntarget)); + + +end \ No newline at end of file diff --git a/Functions/01. Direction finding/df_FFT.m b/Functions/01. Direction finding/df_FFT.m new file mode 100644 index 0000000..42b7d4e --- /dev/null +++ b/Functions/01. Direction finding/df_FFT.m @@ -0,0 +1,40 @@ +function [esti_ang_deg, P_FT, fft_ang_deg] = df_FFT(snapshot, NFFT, fft_ang, array_struct, Ntarget, plot_flag) + +%% FFT algorithm +Nsnap = size(snapshot, 2); + +if isempty(fft_ang) + fft_freq = -180 : 360/NFFT : (180 - 360/NFFT); + fft_ang_deg = asind(fft_freq / 2 * array_struct.lambda_c / array_struct.unit / 180); +end + +% FFT spectrum +zero_padded_input = complex(zeros(NFFT, Nsnap)); +zero_padded_input(array_struct.eff_ch_loc+1, :) = snapshot; + +if Nsnap == 1 + P_FT = abs(flipud(fftshift((1/Nsnap) * ((fft(zero_padded_input, NFFT)))').')).^2; +else + P_FT = abs(flipud(fftshift((1/Nsnap) * (sum(fft(zero_padded_input, NFFT),2))').')).^2; +end + +% Plot +if plot_flag == 1 + figure() + plot(fft_ang_deg, pow2db(abs(P_FT))); + grid on + xlabel('Angle [deg]') + ylabel('Magnitude [dB]') +end + +[pks, pks_ang, ~, ~] = findpeaks(pow2db(abs(P_FT)), fft_ang_deg); +[~, order] = sort(pks, 'descend'); +if isempty(pks) + esti_ang_deg = nan; +else + esti_ang_deg = pks_ang(order(1:Ntarget)); +end + + +end + diff --git a/Functions/01. Direction finding/df_FFT2d_mar510.m b/Functions/01. Direction finding/df_FFT2d_mar510.m new file mode 100644 index 0000000..7ab1e30 --- /dev/null +++ b/Functions/01. Direction finding/df_FFT2d_mar510.m @@ -0,0 +1,62 @@ +function [esti_ang_deg, P_FT, fft_ang_elev_deg, fft_ang_azi_deg] = df_FFT_mar510(zp_snapshot, NFFT_elev, NFFT_azi, fft_ang_elev_deg, fft_ang_azi_deg, array_struct, Ntarget) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Conventioanl Beamforming algorithm for df +% Start : 23.08.25 +% End : 23.08.25 +% developed by Kwanggoo Yeo +% +% Description +% - Input +% - 1) snapshot [matrix], [-] : raw data(snapshot) for df +% - 2) test_ang_elev_deg [vector], [deg] : elevation angle grid for beamforming +% - 3) test_ang_azi_deg [vector], [deg] : Azimuth angle grid for beamforming +% - 4) lambda_c [scalar], [m] : wavelength of center frequency +% - 4) array_struct [struct], [-] : Structure containing virtual array information +% - 5) Ntarget [scalar], [-] : The number of targets in snapshot +% +% - Output +% - 1) esti_ang_deg [matrix], [deg] : The estimated angles by CBF +% - 2) P_CBF [matrix], [linear] : power spectrum of CBF +% +% History +% (23.08.25) Completed +% +% Referece +% - +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +temp_Nsnap = size(zp_snapshot); +Nsnap = temp_Nsnap(end); + +if isempty(fft_ang_elev_deg) + fft_elev_freq = -180 : 360/NFFT_elev : (180 - 360/NFFT_elev); + fft_ang_elev_deg = asind(fft_elev_freq / 2 * array_struct.lambda_c / array_struct.d_unit_elev / 180); +end + +if isempty(fft_ang_azi_deg) + fft_azi_freq = -180 : 360/NFFT_azi : (180 - 360/NFFT_azi); + fft_ang_azi_deg = asind(fft_azi_freq / 2 * array_struct.lambda_c / array_struct.d_unit_azi / 180); +end + +% FFT spectrum +if Nsnap == 1 + P_FT = abs((((1/Nsnap) * ((fft(zero_padded_input, NFFT_azi)))').')).^2; +else + P_FT = abs((fftshift((1/Nsnap) * (sum(fft(zero_padded_input, NFFT_azi),2))').')).^2; +end + + +[pks, pks_ang, ~, ~] = findpeaks(pow2db(abs(P_FT_2D)), fft_ang_deg); +[~, order] = sort(pks, 'descend'); +if isempty(pks) + esti_ang_deg = nan; +else + esti_ang_deg = pks_ang(order(1:Ntarget)); +end + + +end + diff --git a/Functions/01. Direction finding/df_FFT_mar510.m b/Functions/01. Direction finding/df_FFT_mar510.m new file mode 100644 index 0000000..b67d8e7 --- /dev/null +++ b/Functions/01. Direction finding/df_FFT_mar510.m @@ -0,0 +1,60 @@ +function [esti_ang_deg, P_FT, fft_ang_deg] = df_FFT_mar510(snapshot, NFFT, lambda_c, fft_ang, array_struct, Ntarget) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Conventioanl Beamforming algorithm for df +% Start : 23.08.25 +% End : 23.08.25 +% developed by Kwanggoo Yeo +% +% Description +% - Input +% - 1) snapshot [matrix], [-] : raw data(snapshot) for df +% - 2) test_ang_elev_deg [vector], [deg] : elevation angle grid for beamforming +% - 3) test_ang_azi_deg [vector], [deg] : Azimuth angle grid for beamforming +% - 4) lambda_c [scalar], [m] : wavelength of center frequency +% - 4) array_struct [struct], [-] : Structure containing virtual array information +% - 5) Ntarget [scalar], [-] : The number of targets in snapshot +% +% - Output +% - 1) esti_ang_deg [matrix], [deg] : The estimated angles by CBF +% - 2) P_CBF [matrix], [linear] : power spectrum of CBF +% +% History +% (23.08.25) Completed +% +% Referece +% - +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +Nsnap = size(snapshot, 2); + +if isempty(fft_ang) + fft_freq = -180 : 360/NFFT : (180 - 360/NFFT); + fft_ang_deg = asind(fft_freq / 2 * lambda_c / array_struct.d_unit / 180); +end + +% FFT spectrum +zero_padded_input = complex(zeros(NFFT, Nsnap)); +zero_padded_input(flip(abs(array_struct.azi_eff_ch_loc))+1, :) = snapshot; + +if Nsnap == 1 + %P_FT = abs((fftshift((1/Nsnap) * ((fft(zero_padded_input, NFFT)))').')).^2; + P_FT = abs((((1/Nsnap) * ((fft(zero_padded_input, NFFT)))').')).^2; +else + P_FT = abs((fftshift((1/Nsnap) * (sum(fft(zero_padded_input, NFFT),2))').')).^2; +end + + +[pks, pks_ang, ~, ~] = findpeaks(pow2db(abs(P_FT)), fft_ang_deg); +[~, order] = sort(pks, 'descend'); +if isempty(pks) + esti_ang_deg = nan; +else + esti_ang_deg = pks_ang(order(1:Ntarget)); +end + + +end + diff --git a/Functions/01. Direction finding/df_OMP.m b/Functions/01. Direction finding/df_OMP.m new file mode 100644 index 0000000..d44c5ef --- /dev/null +++ b/Functions/01. Direction finding/df_OMP.m @@ -0,0 +1,61 @@ +function [esti_ang_deg, P_OMP, test_ang_deg] = df_OMP(snapshot, fov_deg, unit_ang_deg, array_struct, sparsity_val, plot_flag) + +test_ang_deg = -fov_deg : unit_ang_deg : fov_deg - unit_ang_deg; + +Nant = length(array_struct.eff_ch_loc); + +% Spectrum +dic_mat = zeros(Nant, length(test_ang_deg)); +for ang_idx = 1 : length(test_ang_deg) + sv = exp(-1i * 2 * pi / array_struct.lambda_c * array_struct.eff_ch_loc.' * array_struct.unit * sind(test_ang_deg(ang_idx))); + dic_mat(:, ang_idx) = sv; +end + +mc_val = cal_mutual_coherence(dic_mat, 0) + +[N, K] = size(dic_mat); % N:dim of signal, K: # atoms in dictionary +if (N ~= size(snapshot)) + error('Dimension not matched'); +end + +%% Initializing +x = zeros(K,1); % coefficient (output) +r = snapshot; % residual of b +omega = zeros(sparsity_val,1); % selected support +A_omega = []; % corresponding columns of A +cnt = 0; + +%% Iteration +while (cnt < sparsity_val) % choose sparsity_val atoms + cnt = cnt+1; + x_tmp = zeros(K,1); + inds = setdiff([1:K],omega); % iterate all columns except for the chosen ones + for idx = inds + x_tmp(idx) = dic_mat(:,idx)' * r / norm(dic_mat(:,idx)); % sol of min ||a'x-b|| + end + [~, ichosen] = max(abs(x_tmp)); % choose the maximum + omega(cnt) = ichosen; + A_omega = [A_omega dic_mat(:,ichosen)]; + x_ls = A_omega \ snapshot; % Aomega * x_ls = b + r = snapshot - A_omega * x_ls; % update r +end + +for idx = 1 : sparsity_val + x(omega(idx)) = x_ls(idx); %x_sparse(i).value; +end + +P_OMP = zeros(1, K); +P_OMP(omega) = abs(x_ls).^2; + +% Plot +if plot_flag == 1 + figure() + plot(test_ang_deg, (P_OMP/max(P_OMP))); + grid on + xlabel('Angle [deg]') + ylabel('Normalized Magnitude [linear]') +end + +esti_ang_deg = test_ang_deg(omega); + +end \ No newline at end of file diff --git a/Functions/01. Direction finding/df_Phase_monopulse.m b/Functions/01. Direction finding/df_Phase_monopulse.m new file mode 100644 index 0000000..b9f816f --- /dev/null +++ b/Functions/01. Direction finding/df_Phase_monopulse.m @@ -0,0 +1,5 @@ +function [esti_ang_deg] = df_Phase_monopulse(ch1_phase, ch2_phase, ch_distance) + +esti_ang_deg = asind((ch1_phase - ch2_phase) / 2 / pi / ch_distance); + +end \ No newline at end of file diff --git a/Functions/01. Direction finding/df_cbf.m b/Functions/01. Direction finding/df_cbf.m new file mode 100644 index 0000000..1e7bd8e --- /dev/null +++ b/Functions/01. Direction finding/df_cbf.m @@ -0,0 +1,67 @@ +function [esti_ang_deg, P_CBF] = df_cbf(snapshot, test_ang_elev_deg, test_ang_azi_deg, lambda_c, array_struct, Ntarget) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Conventioanl Beamforming algorithm for df +% Start : 23.08.25 +% End : 23.08.25 +% developed by Kwanggoo Yeo +% +% Description +% - Input +% - 1) snapshot [matrix], [-] : raw data(snapshot) for df +% - 2) test_ang_elev_deg [vector], [deg] : elevation angle grid for beamforming +% - 3) test_ang_azi_deg [vector], [deg] : Azimuth angle grid for beamforming +% - 4) lambda_c [scalar], [m] : wavelength of center frequency +% - 4) array_struct [struct], [-] : Structure containing virtual array information +% - 5) Ntarget [scalar], [-] : The number of targets in snapshot +% +% - Output +% - 1) esti_ang_deg [matrix], [deg] : The estimated angles by CBF +% - 2) P_CBF [matrix], [linear] : power spectrum of CBF +% +% History +% (23.08.25) Completed +% +% Referece +% - +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +% Spectrum +P_CBF = zeros(length(test_ang_elev_deg), length(test_ang_azi_deg)); +for ang_elev_idx = 1 : length(test_ang_elev_deg) + for ang_azi_idx = 1 : length(test_ang_azi_deg) + sv = exp(-1i * 2 * pi / lambda_c * ( array_struct.azi_eff_ch_loc.' * array_struct.lambda_c * cosd(test_ang_elev_deg(ang_elev_idx)) * sind(test_ang_azi_deg(ang_azi_idx)) + array_struct.elev_eff_ch_loc.' * array_struct.lambda_c * sind(test_ang_elev_deg(ang_elev_idx)))); + P_CBF(ang_elev_idx, ang_azi_idx) = abs((sum(sv'*snapshot))).^2; + end +end + +if size(P_CBF, 1) == 1 + [pks, pks_ang, ~, ~] = findpeaks(pow2db(abs(P_CBF)), test_ang_azi_deg); + [~, order] = sort(pks, 'descend'); + if isempty(pks) + esti_ang_deg = nan; + else + esti_ang_deg = pks_ang(order(1:Ntarget)); + end +elseif size(P_CBF, 2) == 1 + [pks, pks_ang, ~, ~] = findpeaks(pow2db(abs(P_CBF)), test_ang_elev_deg); + [~, order] = sort(pks, 'descend'); + if isempty(pks) + esti_ang_deg = nan; + else + esti_ang_deg = pks_ang(order(1:Ntarget)); + end +else + [pks, locs_x, locs_y] = peaks2(pow2db(abs(P_CBF))); + [~, order] = sort(pks, 'descend'); + if isempty(pks) + esti_ang_deg = nan; + else + esti_ang_deg = [test_ang_elev_deg(locs_x(order(1:Ntarget))).' test_ang_azi_deg(locs_y(order(1:Ntarget))).']; + end +end + +end + diff --git a/Functions/01. Direction finding/gen_sig_for_df.m b/Functions/01. Direction finding/gen_sig_for_df.m new file mode 100644 index 0000000..1c7cef3 --- /dev/null +++ b/Functions/01. Direction finding/gen_sig_for_df.m @@ -0,0 +1,124 @@ +function [snapshot, total_SNR_dB, R_nf, r_m_array, sv_mat, noise] = gen_sig_for_df(fc, target_rng, elev_deg, azi_deg, SNR_dB, txarray_loc, rxarray_loc, Nsnap, ch_error, noise_flag) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Generating raw data(snapshot) for direction finding +% Start : 23.06.01 +% End : 23.08.25 +% developed by Kwanggoo Yeo +% +% Description +% - Input +% - 1) fc [vector], [Hz] : Center freqeuncy of radar Tx signal +% - 2) target_rng [vector], [m] : Distance(range) from the radar to each target +% - 3) elev_deg [vector], [deg] : Elevation angle for each taret +% - 4) azi_deg [vector], [deg] : Azimuth angle for each taret +% - 5) SNR_dB [vector], [dB] : SNR for each taret +% - 6) txarray_loc [matrix], [m] : Locations of txarray elements in cartesian coordinate +% - 7) rxarray_loc [matrix], [m] : Locations of rxarray elements in cartesian coordinate +% - 8) Nsnap [scalar], [-] : The number of snapshots to generate +% - 9) ch_error [vector], [-] : gain & phase error for each element in virtual array +% - 10) noise_flag [scalar], [-] : Flag for noise which is included in snapshot (1) or not (0) +% +% - Output +% - 1) snapshot [matrix], [-] : Generated raw data(snapshot) +% - 2) total_SNR_dB [scalar], [dB] : SNR of snapshot, not SNR of signal in single element +% - 3) R_nf [scalar], [m] : Distance of near-field (Frensel Region) +% - 4) r_m_array [matrix], [m] : Distance from the targets to each element in virtual array +% - 5) sv_amt [matrix], [-] : Steering matrix +% - 6) noise [matrix], [-] : Generated noise in each element in virtual array +% +% History +% (23.08.25) Completed +% +% Referece +% - +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% Assumptions +% +% 1. Non-dispersion medium +% 2. Narrowband +% 3. No assumption on Near/Far-field +% 4. Geometry +% 4-1. Cartesian(x-y-z) 좌표계이며, 레이더 시스템의 local coordinate +% 4-2. Elevation angle : x-y 평면 기준으로 위 : + , 아래 : - +% 4-3. Azimuth angle : x-y 평면 상에서, x축의 왼쪽 : +, 오른쪽 : - +% 4-4. Virtual array (tx1, rx1) 채널의 위치를 원점으로 기준(coordinate origin) +% 4-5. x축 : 레이더 boresignt 방향 +% 5. Monostatic radar system + +%% Basic constants settings + +c0 = physconst('LightSpeed'); +lambda_c = c0/fc; +k_c = 2*pi / lambda_c; + +Np = length(azi_deg); % Number of targets + +% ch_error must be column vector +if size(ch_error, 1) == 1 + ch_error = ch_error.'; +end + +%% Array settings + +Nt = size(txarray_loc, 2); +Nr = size(rxarray_loc, 2); + +% Near field region calcuation +array_dist = zeros(Nt*Nr, 1); +for txidx = 1 : Nt + for rxidx = 1 : Nr + array_dist(Nr*(txidx-1) + rxidx) = sqrt(sum((txarray_loc(:,txidx) - rxarray_loc(:,rxidx)).^2)); + end +end +max_len_array = max(array_dist); +R_nf = 2 * max_len_array^2 / lambda_c; + +% Target location calcuation in Cartesian +target_loc = zeros(3, Np); +for tidx = 1 : Np + target_loc(:, tidx) = [target_rng(tidx)*cosd(elev_deg(tidx))*sind(azi_deg(tidx)) ; target_rng(tidx)*cosd(elev_deg(tidx))*cosd(azi_deg(tidx)); target_rng(tidx)*sind(elev_deg(tidx))]; +end + +% Generating steering vector(sv) matrix +r_m_array = complex(zeros(Nt*Nr, Np)); +sv_mat = complex(zeros(Nt*Nr, Np)); + +for target_idx = 1 : Np + for txarray_idx = 1 : Nt + for rxarray_idx = 1 : Nr + r_m_array(Nr*(txarray_idx-1)+rxarray_idx, target_idx) = sqrt(sum((txarray_loc(:, txarray_idx) - target_loc(:, target_idx)).^2)) + sqrt(sum((rxarray_loc(:, rxarray_idx) - target_loc(:, target_idx)).^2)); + sv_mat(Nr*(txarray_idx-1)+rxarray_idx, target_idx) = exp(1i * k_c * r_m_array(Nr*(txarray_idx-1)+rxarray_idx, target_idx)); + end + end + sv_mat(:, target_idx) = sv_mat(:, target_idx) .* ch_error; + sv_mat(:, target_idx) = sv_mat(:, target_idx) * conj(sv_mat(1, target_idx)); +end + + +%% Signal Generation +% Calculating signal complex gain baed on SNR_dB +noise = (sqrt(0.5) * (randn(Nt*Nr, Nsnap) + 1i * randn(Nt*Nr, Nsnap))); +Pn = sum(diag(noise * noise')/Nsnap) / (Nt*Nr); + +sig = zeros(Np, Nsnap); +for target_idx = 1 : Np + Ps = 10^(SNR_dB(target_idx)/10) * Pn; + sig(target_idx, :) = sqrt(Ps) * exp(1i * 2 * pi * rand(1, Nsnap)); +end + +% Signal model +if noise_flag == 1 + snapshot = sv_mat * sig + noise; +else + snapshot = sv_mat * sig; +end + +% Total SNR_dB +Ps_total = sum(diag((sv_mat * sig) * (sv_mat * sig)')/Nsnap) / (Nt*Nr); +total_SNR_dB = 10 * log10( Ps_total / Pn ); + +end + diff --git a/Functions/01. Direction finding/gen_virarray.m b/Functions/01. Direction finding/gen_virarray.m new file mode 100644 index 0000000..a55e407 --- /dev/null +++ b/Functions/01. Direction finding/gen_virarray.m @@ -0,0 +1,36 @@ +function [mimoarray] = gen_virarray(txarray_loc, rxarray_loc) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Generating MIMO virtual array +% Start : 23.08.25 +% End : 23.08.25 +% developed by Kwanggoo Yeo +% +% Description +% - Input +% - 1) txarray_loc [matrix], [m] : Locations of txarray elements in cartesian coordinate +% - 2) rxarray_loc [matrix], [m] : Locations of rxarray elements in cartesian coordinate +% +% - Output +% - 1) mimoarray [matrix], [m] : virtual array +% +% History +% (23.08.25) Completed +% +% Referece +% - +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +Nt = size(txarray_loc,2); +Nr = size(rxarray_loc,2); + +mimoarray = zeros(3, Nt*Nr); +for tx_idx = 1 : Nt + for rx_idx = 1 : Nr + mimoarray(:, Nr*(tx_idx-1) + rx_idx) = txarray_loc(:, tx_idx) + rxarray_loc(:, rx_idx); + end +end + +end + diff --git a/Functions/03. FMCW radar design/SNR_cal.m b/Functions/03. FMCW radar design/SNR_cal.m new file mode 100644 index 0000000..827a0b7 --- /dev/null +++ b/Functions/03. FMCW radar design/SNR_cal.m @@ -0,0 +1,68 @@ +clear all +close all +clc + +c0 = physconst('LightSpeed'); +fc = 76.5e9; +lambda_c = c0 / fc; +k_c = 2 * pi / lambda_c; + +num_tgt = 2; + +Ntar_sample = 300; +tar_loc_x = linspace(-2.5, 2.5, num_tgt); +tar_loc_y = linspace(18, 350, Ntar_sample); +tar_loc_z = zeros(1, num_tgt); + +for tidx = 1 : num_tgt + tar_xyz(:,:,tidx) = [tar_loc_x(tidx)*ones(Ntar_sample,1) tar_loc_y' tar_loc_z(tidx)*ones(Ntar_sample,1)]; + [azimuth, elevation, r] = cart2sph(tar_xyz(:,1,tidx), tar_xyz(:,2,tidx), tar_xyz(:,3,tidx)); + tar_scs(:,:,tidx) = [rad2deg(pi/2 - azimuth) rad2deg(elevation) r]; +end + +Pt = 13; % [dBm] +Pt = Pt - 30; % [dBW] +NF = 12; % [dB] +T0 = 300; % [K] +kb = physconst('Boltzmann'); +BW_int = 40e6; % [Hz] instantaneous Bandwidth + +% Target RCS +rcs = 10; % [dBsm] + +% Noise & Quantization error power +%[~, Qnf] = quanttemp(T0, 12, 'DynamicRange', 52); % ADC bit : 12 [bits], Dynamic Range : 52 [dB] (NXP chip) +%Qnf = 0; +N_bits = 12; +Qnf = -1 * (6.02*N_bits + 10*log10(BW_int) + 1.76) - 30; +Pnq = pow2db(kb * T0) + NF + Qnf + pow2db(BW_int); % [dBW] + +% Antenna patter loading +load('azi_ant_pat.mat'); +% azi_ant_pat(:,2) = azi_ant_pat(:,2) -16.3 + 15; +load('elev_ant_pat.mat'); +% elev_ant_pat(:,2) = elev_ant_pat(:,2) -16.3 + 15; + +for tidx = 1 : num_tgt + tar_ant_gain(:, tidx) = interp1(azi_ant_pat(:,1), azi_ant_pat(:,2), tar_scs(:,1,tidx)); + tar_elev_ant_gain_reduction = interp1(elev_ant_pat(:,1), elev_ant_pat(:,2), 0) - interp1(elev_ant_pat(:,1), elev_ant_pat(:,2), tar_scs(:,2,tidx)); + tar_ant_gain(:, tidx) = tar_ant_gain(:,tidx) - tar_elev_ant_gain_reduction; +end + +% Losses & SP gains +L_sf = 3; % [dB] secondary surface loss +L_ant = 2.5; % [dB] Feeder and Radome loss +L_win = 2.38 + 1.36; % [dB] loss by windowing in 2D +L_straddle = 2.88; % [dB] straddle loss(worst) +L_q = 1; % [dB] Quantization loss +L_sp = 1; % [dB] Other signal processing losses +L_total = L_sf + L_ant + L_win + L_straddle + L_q + L_sp; + +NFFT_R = 512; +NFFT_D = 256; +G_sp = pow2db(NFFT_R * NFFT_D); % [dB] Signal processing gain by 2D-FFT + +% SNR calculation +for tidx = 1 : num_tgt + SNR_set(:,tidx) = Pt + pow2db((lambda_c^2) / (4*pi)^3) + pow2db(1./tar_scs(:,3,tidx).^4) + tar_ant_gain(:,tidx) + tar_ant_gain(:,tidx) + rcs - L_total + G_sp - Pnq; +end diff --git a/Functions/03. FMCW radar design/azi_ant_pat.mat b/Functions/03. FMCW radar design/azi_ant_pat.mat new file mode 100644 index 0000000000000000000000000000000000000000..c1f7fc79ef13b8b49f59ee70540853e2b3b031ac GIT binary patch literal 1003 zcmV_4ZaN@Tb!8w*VQwHXHXt!IIx{pnF)<)AFfueCARr(hARr(hARr(hARr(hARr(h zARr(hARr(hARr(hARr(hARr(B00000000000ZB~{0001E0{{ScoGnrZOp{R%{zWRN zSSkup6ey@2*wViNL`3aJ1Odw`g<=pb;6U6WVnP68#29gaaXVv2cr=Fgd zWlmp}V_uV>O^wi=jrdlJR^w?+7_B- zF^~!Lx~#}RBy)qCZ^#T}`-*pekM+co8N(hOTr)qK(CqGx2tA>%-?}bEk$XRluMdBW zBm!)=X9P(txolC?hY`8fv3~rMP_p2Bl#|`|>BMC-xm9~Oh}fIN>eKyNvQ_}(eTe|~ z()!5@uquba!9AF0eu8&&&msRc@=-7;O1yIOi!E4^~$!1JmIv zLtS{wkNyB2@jJ^_xg6v{S<^~dME2y3Zg@8Ztu^kN=Zk_M$jdDZL2L$9nll|t{}zC< zJMiT!Z1I`1Y;hgosJbd`%>ES#rY9r{0{mu957v)%Js8c_XfWL=arjQ#Em;~5igx6> z(Ez4bl3)OfaO8cY-$1$-2E6;N z>!lOnT;*TaoGb&Izg31C4l^o~2IP3ZY_heiKmj4G?urs8W+u-&)#eLUN20>9sD#X8z`j# z8q{Vs^{#W(V6i~N!4GWz`qX&#Xw$mc0cw0A!4u|;Q)8#1Z$vAvf`9Sb{iabW825Yc zT7TacZ2wD@|Nj+e*VBQ<U=nT_A(72F(f?vbfco>_wAcR^FDh|M5Yg$JETQC~wd% z)tpd!p}4Ws3(X|ro2DL@EXWF2cq7ff-Ouy!U ZZfLRB!p|vh+PclTah`70wtqUHlF)n<&v5_% literal 0 HcmV?d00001 diff --git a/Functions/03. FMCW radar design/beat2rng.m b/Functions/03. FMCW radar design/beat2rng.m new file mode 100644 index 0000000..eb34e90 --- /dev/null +++ b/Functions/03. FMCW radar design/beat2rng.m @@ -0,0 +1,11 @@ +function [rng] = beat2rng(f_beat, f_slope, c0) + +% Objective : Converting beat frequency [Hz] to range [m] + +if isempty(c0) + c0 = physconst('LightSpeed'); +end + +rng = c0 * f_beat / f_slope / 2; + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/bw2rngsep.m b/Functions/03. FMCW radar design/bw2rngsep.m new file mode 100644 index 0000000..f6a61a9 --- /dev/null +++ b/Functions/03. FMCW radar design/bw2rngsep.m @@ -0,0 +1,11 @@ +function [rng_sep] = bw2rngsep(TxBw, rb, c0) + +% Objective : Converting range resolution to required 3-dB bandwidth needed to distinguish two targets separated by the range specified in r [m] + +if isempty(c0) + c0 = physconst('LightSpeed'); +end + +rng_sep = (c0 * rb) ./ (2*TxBw); + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/cal_nf.m b/Functions/03. FMCW radar design/cal_nf.m new file mode 100644 index 0000000..d32a1a2 --- /dev/null +++ b/Functions/03. FMCW radar design/cal_nf.m @@ -0,0 +1,38 @@ +function [total_nf, total_gain, total_temp] = cal_nf(nf_set, gain_set, reftemp) + +% Check inputs +narginchk(2,3); + +% Validate +validateattributes(nf_set,{'double'}, {'nonnan','nonempty','real', ... + 'vector','nonnegative'}, 'noisefigure', 'NF'); +numNF = numel(nf_set); +validateattributes(gain_set,{'double'}, {'nonnan','nonempty','real', ... + 'vector','numel',numNF}, 'noisefigure', 'G'); +NFcol = db2pow(nf_set(:)); +Gcol = db2pow(gain_set(:)); + +% Check for temperature input +if nargin < 3 + reftemp = 290; % Standard noise temperature (Kelvin) +else + validateattributes(reftemp,{'double'}, {'finite','nonempty','real', ... + 'nonnegative','scalar'}, 'noisefigure', 'REFTEMP'); +end + +% Calculate total gain +total_gain = sum(gain_set(:),1); % dB + +% Calculate cascaded noise figure +if numel(NFcol) > 1 + cnfLinear = NFcol(1) + sum((NFcol(2:end) - 1)./cumprod(Gcol(1:end-1),1),1); % Linear + total_nf = 10*log10(cnfLinear); % dB +else + cnfLinear = NFcol(1); + total_nf = nf(1); +end + +% Calculate cascaded noise temperature +total_temp = reftemp*(cnfLinear); % Kelvin + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/cal_rdcoupling.m b/Functions/03. FMCW radar design/cal_rdcoupling.m new file mode 100644 index 0000000..aaddcc3 --- /dev/null +++ b/Functions/03. FMCW radar design/cal_rdcoupling.m @@ -0,0 +1,11 @@ +function [deltaR] = cal_rdcoupling(dop, f_slope, c0) + +% Objective : Calculating range offset [m] due to Doppler shift [Hz] in a LFM signal + +if isempty(c0) + c0 = physconst('LightSpeed'); +end + +deltaR = -c0 * dop / (2 * f_slope); + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/cal_rtmcoupling.m b/Functions/03. FMCW radar design/cal_rtmcoupling.m new file mode 100644 index 0000000..a7349fa --- /dev/null +++ b/Functions/03. FMCW radar design/cal_rtmcoupling.m @@ -0,0 +1,7 @@ +function [deltaR] = cal_rtmcoupling(spd, PRI, m) + +% Objective : Calculating range offset [m] due to Target Motion during sweeps in a LFM signal + +deltaR = -m * PRI * spd; + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/dechirping.m b/Functions/03. FMCW radar design/dechirping.m new file mode 100644 index 0000000..52a9a6a --- /dev/null +++ b/Functions/03. FMCW radar design/dechirping.m @@ -0,0 +1,10 @@ +function [sig_out] = dechirping(sig_in, sig_ref) + +% Objecrtive : Mixing the incoming signal, sig_in, with the reference +% signal, sig_ref. + +sig_ref = cast(sig_ref, class(sig_in)); + +sig_out = bsxfun(@times, conj(sig_ref), sig_in); + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/dop2spd.m b/Functions/03. FMCW radar design/dop2spd.m new file mode 100644 index 0000000..84271c7 --- /dev/null +++ b/Functions/03. FMCW radar design/dop2spd.m @@ -0,0 +1,7 @@ +function [spd] = dop2spd(dop, lambda) + +% Objective : Converting Doppler shift [Hz] to Radial Speed [m/s] for one-way propagation + +spd = dop * lambda; + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/elev_ant_pat.mat b/Functions/03. FMCW radar design/elev_ant_pat.mat new file mode 100644 index 0000000000000000000000000000000000000000..fb48ae1a6c02410d374112497d771acaa107504c GIT binary patch literal 1064 zcmV+@1lRjbK~zjZLLfCRFd$7qR4ry{Y-KDUP;6mzW^ZzBIv`L(S4mDbG%O%Pa%Ew3 zWn>_4ZaN@Tb!8w*VQwHXHXt!IIy5pmGczDEFfueCARr(hARr(hARr(hARr(hARr(h zARr(hARr(hARr(hARr(hARr(B00000000000ZB~{0001>0{{ScoIR8aNK|PQ$8Tk> znWU)Xl+!d#*Ll@AGdiPVsdMg_iKRBE>sO;Ws^f_5F`x z^%kp_STSPFs+_f;tb9&&Syg#8+4U3azX(NnR?=pwNgj*G-t=5sS`W-5x>QjJl?MlVNvF4ZF)@zL(?B!wtUubSHAFYRKuIQt8$`~ zG>lxiF>mb`V7u#pDfFQ~U`+23zYbg+io2){bHcneDYIlj+C`nfsYgE zK!1oP;npO8g90=eAgSxphlvl2FepNJer<-o&o|~Dn`Z%%^U(wq5?c!HFPWt3c2?lJ zZLs3!6V)4TIbfbBf7yzW@K%)nWM3+9-X37S9-y-BMWg}ikPb54G9PU3G#{G#zaNH? zfW-kATJyA|?VnK~`7FN;IJGhE5HbPY8nkX~zpIL$Gg~ zaNH3BFA38bU_H$lI7GPl_n=?r>0>=Fq_d8ghD)i?DPEnUSJ&N*zvPf&F6MF9fh0%b z3p+~v>~=gy7`s*G8?7?O`r%0JdM7IdNyf_-tKQc(?u#vUlpLZRD7>{4G&=C&1$vy} zcrPjZIX4P*gZ{4`uA{-j{XO1=M7gBR#&c9V45_;sPnq<+0VSWh1SBg{O|0He-`XXI zmEVF*!u@J37RJ$_Y>H<-+HmH7LB)b zMg#npF|6y1IMF%3jX;`?{gJ79y-xievGLr!OTo33;=^HUn9rUn5?tju)vPG{^`SlN iM{;$hAzde3E$Pg?$%i&-uq7F|m@bBO0sIa(@kHn3V)kqR literal 0 HcmV?d00001 diff --git a/Functions/03. FMCW radar design/gen_fastramp_fmcw.m b/Functions/03. FMCW radar design/gen_fastramp_fmcw.m new file mode 100644 index 0000000..3dc219c --- /dev/null +++ b/Functions/03. FMCW radar design/gen_fastramp_fmcw.m @@ -0,0 +1,68 @@ +function [waveform_seq] = gen_fastramp_fmcw(timing_struct, fs, f0, TxBw, NumChirps, prop_delay, phase_initial, opt) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Generating Fast Ramp Linear FMCW waveform +% Start : 23.12.22 +% End : xx.xx.xx +% developed by Kwanggoo Yeo +% +% Description +% - Input +% - 1) timing_struct [struct], [-] : Timing information for one chirp +% - 1-1) T_dwell [scalar], [sec] : Time length for idle +% - 1-2) T_settle [scalar], [sec] : Time length for ramp start but not acquisition +% - 1-3) T_jumpback [scalar], [sec] : Time length for go back to the start frequency +% - 1-4) T_reset [scalar], [sec] : Time length for reset +% - 1-5) T_acq [scalar], [sec] : Time length for acqusition +% - 2) fs [scalar], [Hz] : Sampling rate +% - 3) f0 [scalar], [Hz] : Start frequency +% - 4) TxBw [scalar], [Hz] : Waveform Bandwidth +% - 5) NumChirps [scalar], [-] : The number of chirps in one frame +% - 6) prop_delay [scalar], [sec] : Time delay by range between radar and target +% - 7) phase_ini [scalar], [deg] : Initial phase of Tx waveform in deg +% +% - Output +% - 1) waveform_seq [struct], [-] : Struct containing waveform information +% - 1-1) T_chirp [scalar], [sec] : Time length for one chirp +% - 1-2) T_frame [scalar], [sec] : Time length for one frame +% - 1-3) timeline [vec], [sec] : Time index for Fast ramp Linear FMCW sequence +% - 1-4) waveform [matrix], [-] : Fast ramp Linear FMCW sequence +% +% History +% (23.12.22) Start +% +% Referece +% - +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +T_dwell = timing_struct.T_dwell; +T_settle = timing_struct.T_settle; +T_acq = timing_struct.T_acq; +T_jumpback = timing_struct.T_jumpback; +T_reset = timing_struct.T_reset; +T_idle = timing_struct.T_idle; + +waveform_seq.T_chirp = T_dwell + T_settle + T_acq + T_jumpback + T_reset + T_idle; + +waveform_seq.T_frame = waveform_seq.T_chirp * NumChirps; + +waveform_seq.f_slope = TxBw / T_acq; + +waveform_seq.timeline = (0 : 1/fs : (T_acq - 1/fs)) + T_dwell + T_settle + prop_delay; + +if opt == 0 + waveform_seq.waveform = cos(2 * pi * f0 * waveform_seq.timeline + 2 * pi * (waveform_seq.f_slope/2 * waveform_seq.timeline.^2) + phase_initial); + waveform_seq.waveform = waveform_seq.waveform / sqrt(norm(waveform_seq.waveform)^2 / length(waveform_seq.waveform)); +elseif opt == 1 + waveform_seq.waveform = exp(1i * 2 * pi * f0 * waveform_seq.timeline + 1i * 2 * pi * (waveform_seq.f_slope/2 * waveform_seq.timeline.^2) + 1i* phase_initial); + waveform_seq.waveform = waveform_seq.waveform / sqrt(norm(waveform_seq.waveform)^2 / length(waveform_seq.waveform)); +else + disp('Opt error! : Real waveform(opt = 0) or Complex waveform(opt = 1)') + return; +end + + + +end + diff --git a/Functions/03. FMCW radar design/gen_fastramp_fmcw2.m b/Functions/03. FMCW radar design/gen_fastramp_fmcw2.m new file mode 100644 index 0000000..59af687 --- /dev/null +++ b/Functions/03. FMCW radar design/gen_fastramp_fmcw2.m @@ -0,0 +1,71 @@ +function [waveform_seq] = gen_fastramp_fmcw2(timing_struct, fs, f0, TxBw, ChirpIndex, NumChirps, prop_delay, phase_initial, opt) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Generating m-th Fast Ramp Linear FMCW waveform with considering different PRI +% Start : 23.12.22 +% End : xx.xx.xx +% developed by Kwanggoo Yeo +% +% Description +% - Input +% - 1) timing_struct [struct], [-] : Timing information for one chirp +% - 1-1) T_dwell [scalar], [sec] : Time length for idle +% - 1-2) T_settle [scalar], [sec] : Time length for ramp start but not acquisition +% - 1-3) T_jumpback [scalar], [sec] : Time length for go back to the start frequency +% - 1-4) T_reset [scalar], [sec] : Time length for reset +% - 1-5) T_acq [scalar], [sec] : Time length for acqusition +% - 2) fs [scalar], [Hz] : Sampling rate +% - 3) f0 [scalar], [Hz] : Start frequency +% - 4) TxBw [scalar], [Hz] : Waveform Bandwidth +% - 5) ChirpIndex [scalar], [-] : m-th chirp index ( m= 1, 2, 3, ..., M ) +% - 6) NumChirps [scalar], [-] : The number of chirps in one frame +% - 7) prop_delay [scalar], [sec] : Time delay by range between radar and target +% - 8) phase_ini [scalar], [deg] : Initial phase of Tx waveform in deg +% +% - Output +% - 1) waveform_seq [struct], [-] : Struct containing waveform information +% - 1-1) T_chirp [scalar], [sec] : Time length for one chirp +% - 1-2) T_frame [scalar], [sec] : Time length for one frame +% - 1-3) timeline [vec], [sec] : Time index for Fast ramp Linear FMCW sequence +% - 1-4) waveform [matrix], [-] : Fast ramp Linear FMCW sequence +% +% History +% (23.12.22) Start +% +% Referece +% - +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +T_dwell = timing_struct.T_dwell; +T_settle = timing_struct.T_settle; +T_acq = timing_struct.T_acq; +T_jumpback = timing_struct.T_jumpback; +T_reset = timing_struct.T_reset; +T_idle = timing_struct.T_idle; + +waveform_seq.T_chirp = T_dwell + T_settle + T_acq + T_jumpback + T_reset + T_idle; + +waveform_seq.T_frame = waveform_seq.T_chirp * NumChirps; + +waveform_seq.f_slope = TxBw / T_acq; + +waveform_seq.timeline = T_dwell + T_settle + (0 : 1/fs : (T_acq - 1/fs)) + ((ChirpIndex-1) * waveform_seq.T_chirp) + prop_delay; + + + +if opt == 0 + waveform_seq.waveform = cos(2 * pi * f0 * waveform_seq.timeline + 2 * pi * (waveform_seq.f_slope/2 * waveform_seq.timeline.^2) + phase_initial); + waveform_seq.waveform = waveform_seq.waveform / sqrt(norm(waveform_seq.waveform)^2 / length(waveform_seq.waveform)); +elseif opt == 1 + waveform_seq.waveform = exp(1i * 2 * pi * f0 * waveform_seq.timeline + 1i * 2 * pi * (waveform_seq.f_slope/2 * waveform_seq.timeline.^2) + 1i* phase_initial); + waveform_seq.waveform = waveform_seq.waveform / sqrt(norm(waveform_seq.waveform)^2 / length(waveform_seq.waveform)); +else + disp('Opt error! : Real waveform(opt = 0) or Complex waveform(opt = 1)') + return; +end + + + +end + diff --git a/Functions/03. FMCW radar design/rng2beat.m b/Functions/03. FMCW radar design/rng2beat.m new file mode 100644 index 0000000..c2028a7 --- /dev/null +++ b/Functions/03. FMCW radar design/rng2beat.m @@ -0,0 +1,15 @@ +function [f_beat] = rng2beat(rng, f_slope, c0) + +% Objective : Converting the range [m] of a dechirped linear FMCW signal to +% its corresponding range, beat frequency [Hz] +% +% f_beat = (2 * f_slope * rng) / c0 +% + +if isempty(c0) + c0 = physconst('LightSpeed'); +end + +f_beat = 2 * rng / c0 * f_slope; + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/rng2time.m b/Functions/03. FMCW radar design/rng2time.m new file mode 100644 index 0000000..7ae1c5f --- /dev/null +++ b/Functions/03. FMCW radar design/rng2time.m @@ -0,0 +1,11 @@ +function [t] = rng2time(r, c0) + +% Objective : Calculating the time that a signal takes to propagate given range, r + +if isempty(c0) + c0 = physconst('LightSpeed'); +end + +t = r / c0; + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/rngres2bw.m b/Functions/03. FMCW radar design/rngres2bw.m new file mode 100644 index 0000000..151fed2 --- /dev/null +++ b/Functions/03. FMCW radar design/rngres2bw.m @@ -0,0 +1,11 @@ +function [req_bw] = rngres2bw(rng_res, rb, c0) + +% Objective : Converting range resolution to required 3-dB bandwidth needed to distinguish two targets separated by the range specified in r [m] + +if isempty(c0) + c0 = physconst('LightSpeed'); +end + +req_bw = (c0 * rb) ./ (2*rng_res); + +end \ No newline at end of file diff --git a/Functions/03. FMCW radar design/spd2dop.m b/Functions/03. FMCW radar design/spd2dop.m new file mode 100644 index 0000000..6587992 --- /dev/null +++ b/Functions/03. FMCW radar design/spd2dop.m @@ -0,0 +1,7 @@ +function [dop] = spd2dop(spd, lambda) + +% Objective : Converting the speed [m/s] to the corresponding Doppler frequency shift [Hz] for one-way propagation + +dop = spd / lambda; + +end \ No newline at end of file diff --git a/Functions/04. Target modelling/cal_rts_attenuation.m b/Functions/04. Target modelling/cal_rts_attenuation.m new file mode 100644 index 0000000..e59465a --- /dev/null +++ b/Functions/04. Target modelling/cal_rts_attenuation.m @@ -0,0 +1,10 @@ +function [pow_attenuation_dbm, amp_attenuation_dbm] = cal_rts_attenuation(range_btw_ant_and_radar, gain_Tx_dB, gain_Rx_dB, range_sim_m, lambda_c, RCS_sim_dbsm) + +gain_Tx_lin = 10^(gain_Tx_dB/10); +gain_Rx_lin = 10^(gain_Rx_dB/10); +RCS_sim_lin = 10^(RCS_sim_dbsm/10); +pow_attenuation_dbm = 10*log10(RCS_sim_lin/range_sim_m^4 * 4 * pi * range_btw_ant_and_radar^4 / gain_Tx_lin / gain_Rx_lin / lambda_c^2) -30; +amp_attenuation_dbm = 10*log10(sqrt(RCS_sim_lin/range_sim_m^4 * 4 * pi * range_btw_ant_and_radar^4 / gain_Tx_lin / gain_Rx_lin / lambda_c^2)) -30; + +end + diff --git a/Functions/04. Target modelling/gen_tgt.m b/Functions/04. Target modelling/gen_tgt.m new file mode 100644 index 0000000..bf87d89 --- /dev/null +++ b/Functions/04. Target modelling/gen_tgt.m @@ -0,0 +1,10 @@ +function [tgt] = gen_tgt(pos, vel, acc, rcs, str_class) + +tgt.pos = pos; +tgt.vel = vel; +tgt.acc = acc; +tgt.rcs = rcs; +tgt.class = str_class; + +end + diff --git a/System_Analysis_v2_0.m b/System_Analysis_v2_0.m new file mode 100644 index 0000000..e7d7ec8 --- /dev/null +++ b/System_Analysis_v2_0.m @@ -0,0 +1,666 @@ +clear variables; +close all; +clc + +addpath(genpath(fullfile(pwd, 'Functions'))); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% 좌표계 : ego car 기준 앞뒤 : x 축, 좌우 : y 축, 위아래 : z 축 +% +% +% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% Control Panel +opt_waveform = 1; % 0 : cos waveform, 1 : exp waveform +opt_interference = 0; +opt_clutter = 0; + +opt_tgt = 0; % 0 : Point targets, 1 : Realistic targerts +opt_noise = 1; +opt_filter = 1; +opt_window = 1; + +flag_plot_physical_array = 0; +flag_plot_virtual_array = 0; +flag_plot_RVM = 1; +flag_plot_noise_esti = 0; +flag_plot_birdview = 1; + +%% Basic Parameters +c0 = physconst('Lightspeed'); % Light speed [m/s] +kb = physconst('Boltzmann'); % Boltzmann Constant +T0 = 290; % Room Temperature [K] + +%% System Parameters +fc = 76.75e9; % Center frequency [Hz] +lambda_c = c0 / fc; % wavelength for center frequency [m] +fs = 20e6; % ADC sampling rate [Hz] +NumChirps = 256; % Number of Chirps (= # of samples in slow-time) +NumSamples = 1024; % Number of samples in fast-time +rngFFTlen = 2^(nextpow2(NumSamples)); % FFT length in fast-time +velFFTlen = 2^(nextpow2(NumChirps)); % FFT length in slow-time + +maxrngFFTidx = 13/16 * rngFFTlen / 2; % Max Range index in FFT + +%% Array Design +Nt = 3; % Number of Tx +Nr = 4; % Number of Rx + +u_azi = 1.96e-3; % [m] +u_elv = 5.488e-3; % [m] + +txarray_loc = [ 0 0 0 ; + 0 -8*u_azi 0 ; + 0 -15*u_azi u_elv;]'; +rxarray_loc = [ 0 0 10*u_azi; + 0 -5*u_azi 10*u_azi; + 0 -9*u_azi 10*u_azi; + 0 -15*u_azi 10*u_azi;]'; + +[mimoarray] = gen_virarray(txarray_loc, rxarray_loc); + +if flag_plot_physical_array == 1 + ant_width = 2.55e-3; + ant_height = 10.77e-3; + + figure + for eidx = 1 : Nt + scatter((txarray_loc(2, eidx)+ant_width/2)*1e3, (txarray_loc(3, eidx)+ant_height/2)*1e3, 'ro'); + hold on + rectangle('Position', 1e3*[(txarray_loc(2, eidx)), (txarray_loc(3, eidx)), ant_width, ant_height]); + hold on + text(txarray_loc(2,eidx)*1e3, txarray_loc(3,eidx)*1e3, ['#', num2str(eidx)], 'Color', 'red', 'HorizontalAlignment', 'center'); + end + grid on + + for eidx = 1 : Nr + scatter((rxarray_loc(2, eidx)+ant_width/2)*1e3, (rxarray_loc(3, eidx)+ant_height/2)*1e3, 'bo'); + hold on + rectangle('Position', 1e3*[(rxarray_loc(2, eidx)), (rxarray_loc(3, eidx)), ant_width, ant_height]); + hold on + text(rxarray_loc(2,eidx)*1e3, rxarray_loc(3, eidx)*1e3, ['#', num2str(eidx)], 'Color', 'blue', 'HorizontalAlignment', 'center'); + end + xlabel('y-axis [mm]'); + ylabel('y-axis [mm]'); + title('Physical array locations'); +end + +if flag_plot_virtual_array == 1 + figure + scatter(mimoarray(2,:) / lambda_c, mimoarray(3,:) / lambda_c, 'o') + grid on + xlabel('y-axis [\lambda]'); + ylabel('z-axis [\lambda]'); + title('Virtual Array Positions'); +end + +%% Environments +% 00. Radar platform +posR = [0 0 0].'; % radar position +velR = [0 0 0].'; % radar velocity +accR = [0 0 0].'; % radar acceleration + +% 01. Target +tgt_rng = [10]; +tgt_azi = [0]; +tgt_elv = [0]; +tgt_pos = [tgt_rng.*cosd(tgt_azi).*cosd(tgt_elv) tgt_rng.*sind(tgt_azi).*cosd(tgt_elv) tgt_rng.*sind(tgt_elv)].'; + +%tgt_pos = [150 0 0].'; +tgt_vel = [5 0 0].'; +tgt_acc = [0 0 0].'; +tgt_rcs = [20]; +tgt_class = 'TCR'; +tgt = gen_tgt(tgt_pos, tgt_vel, tgt_acc, tgt_rcs, 'TCR'); +tgt.num = size(tgt.rcs, 1); + +for tgt_idx = 1 : tgt.num + [tgt_azi_rad(tgt_idx), tgt_elv_rad(tgt_idx), ~] = cart2sph(tgt_pos(1, tgt_idx), tgt_pos(2, tgt_idx), tgt_pos(3, tgt_idx)); + tgt_azi_deg(tgt_idx, 1) = rad2deg(tgt_azi_rad(tgt_idx)); + tgt_elv_deg(tgt_idx, 1) = rad2deg(tgt_elv_rad(tgt_idx)); + tgt_radi_vel(tgt_idx, 1) = norm((tgt_vel(:, tgt_idx) - velR) .* (tgt_pos(:, tgt_idx) - posR)/norm((tgt_pos(:, tgt_idx) - posR))); +end + +tgt_table_truth = table([tgt_rng tgt_radi_vel tgt_azi_deg tgt_elv_deg]', 'VariableNames', {'Tgt.'}, 'RowNames', {'Rng[m]', 'Vel,[m/s]', 'Azi.[deg]','Elv.[deg]'}); +disp(tgt_table_truth); + +% 02. Clutter + +% 03. Interference + +%% Transmitter +% 00. Sampling interval +ts = 1/fs; + +% 00. Transmitter Parameters +Ptx_dBm = 11; + +% 01. Waveform generation +TxBw = 344e6; +timing_struct.T_dwell = 0e-6; +timing_struct.T_settle = 3e-6; +timing_struct.T_jumpback = 0.5e-6; +timing_struct.T_reset = 1e-6; +timing_struct.T_acq = 51.2e-6; +tx_phase_initial = 0; +timing_struct.T_idle = 9.5e-6; %[1e-6 9.5e-6 18e-6]; + +% timing_struct.T_dwell = 0e-6; +% timing_struct.T_settle = 0e-6; +% timing_struct.T_jumpback = 0e-6; +% timing_struct.T_reset = 0e-6; +% timing_struct.T_acq = 51.2e-6; +% tx_phase_initial = 0; +% timing_struct.T_idle = 0e-6; %[1e-6 9.5e-6 18e-6]; +timing_struct.PRI = timing_struct.T_dwell + timing_struct.T_settle + timing_struct.T_jumpback + timing_struct.T_reset + timing_struct.T_acq + timing_struct.T_idle; +%[56.7e-6 65.2e-6 73.7e-6]; + + +% [24.08.22] Lowpass filter 구현 때문에 2배 oversampling 함. +for m = 1 : NumChirps + %Tx_wave = gen_fastramp_fmcw2(timing_struct, 2*fs, fc-TxBw/2, TxBw, m, NumChirps, 0, tx_phase_initial, opt_waveform); + Tx_wave = gen_fastramp_fmcw(timing_struct, 2*fs, fc-TxBw/2, TxBw, NumChirps, 0, tx_phase_initial, opt_waveform); + Ref_wave(:,m) = Tx_wave.waveform; +end +Ref_wave = repmat(Ref_wave, 1, 1, Nr); + +% 02. DDMA +DDMA_freq = [0 1/4 2/4] * 1/Tx_wave.T_chirp; +DDMA_idx = DDMA_freq * Tx_wave.T_chirp * velFFTlen; + +%% System Parameters (Analysis) +rng_max = beat2rng(fs/2*13/16, Tx_wave.f_slope, c0); +rng_res = beat2rng(fs/rngFFTlen, Tx_wave.f_slope, c0); +rng_sep = bw2rngsep(TxBw, 1, c0); + +vel_max = 0.5*dop2spd(1/(Tx_wave.T_chirp)/2, lambda_c); +vel_res = dop2spd(1/Tx_wave.T_frame/velFFTlen, lambda_c); +vel_sep = dop2spd(1/Tx_wave.T_frame/velFFTlen, lambda_c); + +% figure(1) +% subplot(211); plot(Tx_wave.timeline,real(Tx_wave.waveform)); +% xlabel('Time (s)'); ylabel('Amplitude (v)'); +% title('FMCW signal'); axis tight; +% subplot(212); spectrogram(Tx_wave.waveform,32,16,32,fs,'yaxis'); +% title('FMCW signal spectrogram'); + +%% Antenna (Tx) +load('azi_ant_pat.mat'); +load('elev_ant_pat.mat'); + +azi_ant_pat(:,2) = azi_ant_pat(:,2) - 4; +elev_ant_pat(:,2) = elev_ant_pat(:,2) - 4; + +%% Propagation (From Tx ant to Rx ant) +delayed_Tx_sig = zeros(length(Tx_wave.waveform), NumChirps, Nr); + +for tgtidx = 1 : tgt.num + + for m = 1 : NumChirps + + if m > 1 + velR = velR + accR * timing_struct.PRI; + posR = posR + velR * timing_struct.PRI; + tgt.pos(:,tgtidx) = tgt.pos(:,tgtidx) + tgt.vel(:,tgtidx) * timing_struct.PRI; + tgt.vel(:,tgtidx) = tgt.vel(:,tgtidx) + tgt.acc(:,tgtidx) * timing_struct.PRI; + end + + % RF out + powVar_dB = (Ptx_dBm - 30) * ones(Nt, Nr); % -30 : dBm -> dBW + for txidx = 1 : Nt + tgt_loc_vec = tgt.pos(:, tgtidx) - (posR + txarray_loc(:, txidx)); + [tgt_azi_dod, tgt_elev_dod, tgt_rng_dod] = cart2sph(tgt_loc_vec(1),tgt_loc_vec(2),tgt_loc_vec(3)); + target_prop_time_dod = rng2time(tgt_rng_dod, c0); + tx_ant_gain_azi = interp1(azi_ant_pat(:,1), azi_ant_pat(:,2), rad2deg(tgt_azi_dod)); + tx_ant_gain_elev_reduction = max(elev_ant_pat(:,2)) - interp1(elev_ant_pat(:,1), elev_ant_pat(:,2), rad2deg(tgt_elev_dod)); + tx_ant_gain_db = tx_ant_gain_azi - tx_ant_gain_elev_reduction; + + % Tx ant gain + powVar_dB(txidx,:) = powVar_dB(txidx,:) + tx_ant_gain_db; + + % Freespace loss (Radar to Target) + powVar_dB(txidx,:) = powVar_dB(txidx,:) + pow2db(1/(4*pi*tgt_rng_dod^2)); + + % RCS + powVar_dB(txidx,:) = powVar_dB(txidx,:) + tgt.rcs(tgtidx); + + for rxidx = 1 : Nr + tgt_loc_vec = tgt.pos(:, tgtidx) - (posR + rxarray_loc(:, rxidx)); + [tgt_azi_doa, tgt_elev_doa, tgt_rng_doa] = cart2sph(tgt_loc_vec(1),tgt_loc_vec(2),tgt_loc_vec(3)); + target_prop_time_doa = rng2time(tgt_rng_doa, c0); + rx_ant_gain_azi = interp1(azi_ant_pat(:,1), azi_ant_pat(:,2), rad2deg(tgt_azi_doa)); + rx_ant_gain_elev_reduction = max(elev_ant_pat(:,2)) - interp1(elev_ant_pat(:,1), elev_ant_pat(:,2), rad2deg(tgt_elev_doa)); + rx_ant_gain_db = rx_ant_gain_azi - rx_ant_gain_elev_reduction; + + % Freespace loss (Target to Radar) + powVar_dB(txidx, rxidx) = powVar_dB(txidx, rxidx) + pow2db(1/(4*pi*tgt_rng_doa^2)); + + % Rx ant gain + powVar_dB(txidx, rxidx) = powVar_dB(txidx, rxidx) + rx_ant_gain_db + pow2db(lambda_c^2/(4*pi)); + + target_prop_time = target_prop_time_dod + target_prop_time_doa; + % [24.08.22] Lowpass filter 구현 때문에 2배 oversampling 함. + %delayed_wave = gen_fastramp_fmcw2(timing_struct, 2*fs, fc-TxBw/2, TxBw, m, NumChirps, target_prop_time, tx_phase_initial, opt_waveform); + delayed_wave = gen_fastramp_fmcw(timing_struct, 2*fs, fc-TxBw/2, TxBw, NumChirps, target_prop_time, tx_phase_initial, opt_waveform); + delayed_Tx_sig(:,m,rxidx) = delayed_Tx_sig(:,m,rxidx) + sqrt(db2pow(powVar_dB(txidx, rxidx))) * (delayed_wave.waveform).' .* exp(-1i * 2 * pi * DDMA_freq(txidx) * ((m-1) * delayed_wave.T_chirp)); + end + end + end +end + + +%% Antenna (Rx) +rxsig_wo_n = delayed_Tx_sig; + +%% Receiver +% Compute the cascaded noise figure and total gain of a receiver system. The system has seven stages, with these values: +% 01. LNA with a noise figure of 1.0 dB and a gain of 15.0 dB +% 02. Mixer with a noise figure of 5.0 dB and a gain of –7.0 dB +% 03. HPF1 with a noise figure of 0.5 dB and a gain of –0.5 dB +% 04. IF VGA1 with a noise figure of 0.6 dB and a gain of –15.0 dB +% 05. HPF2 with a noise figure of 0.5 dB and a gain of –0.5 dB +% 06. IF VGA2 with a noise figure of 0.6 dB and a gain of –15.0 dB +% 07. LPF1 with a noise figure of 0.6 dB and a gain of –15.0 dB +% 08. Buffer with a noise figure of 1.0 dB and a gain of –1.0 dB +% 09. LPF2 with a noise figure of 0.6 dB and a gain of –15.0 dB +% 10. ADC with Vpp : 1.2 [V], Bit : 12 [bit] + +% 00. Receiver Parameters +NF_dB_LNA = 11.5; +NF_dB_Mixer = 12; +NF_dB_HPF1 = 0.5; +NF_dB_VGA1 = 20; +NF_dB_HPF2 = 0.5; +NF_dB_VGA2 = 20; +NF_dB_LPF1 = 0.5; +NF_dB_Buffer = 30; +NF_dB_LPF2 = 0.5; + +G_dB_LNA = 20; +G_dB_Mixer = 0; +G_dB_HPF1 = -0.5; +G_dB_VGA1 = 13.5; +G_dB_HPF2 = -0.5; +G_dB_VGA2 = 13.5; +G_dB_LPF1 = -0.5; +G_dB_Buffer = 0; +G_dB_LPF2 = -0.5; + + +nf = [NF_dB_LNA NF_dB_Mixer NF_dB_HPF1 NF_dB_VGA1 NF_dB_HPF2, NF_dB_VGA2 NF_dB_LPF1 NF_dB_Buffer NF_dB_LPF2]; +g = [G_dB_LNA G_dB_Mixer G_dB_HPF1 G_dB_VGA1 G_dB_HPF2, G_dB_VGA2 G_dB_LPF1 G_dB_Buffer G_dB_LPF2]; + + +% nf = [4.0 0.5 5.0 1.0 0.6 1.0 6.0]; +% g = [20.0 -0.5 -7.0 -1.0 21.75 21.75 -5.0]; + +[cnf,ng] = noisefigure(nf,g) +% NF_dB2 = cal_nf(nf, g, T0); +NF_dB = 12; + +% 01. LNA +NF_dB_LNA = 12; +if opt_noise == 0 + Pn = 0; +elseif opt_noise == 1 + Pn = kb * T0 * fs * db2pow(NF_dB_LNA); +end + +thermal_noise_LNA = sqrt(Pn) * randn([size(rxsig_wo_n)]); + +% 02. Mixing +% 2-1. Product +mixed_sig = zeros(size(rxsig_wo_n)); +for rxidx = 1 : Nr + mixed_sig(:,:,rxidx) = dechirping(rxsig_wo_n(:,:,rxidx), Ref_wave(:,:,rxidx)); + mixed_noise(:,:,rxidx) = dechirping(thermal_noise_LNA(:,:,rxidx), Ref_wave(:,:,rxidx)); +end + +% 2-2. High Pass Filter +% Second-order, -6 dB, (100 kHz ~ 3.2 MHz) : NXP (tef82xx) +f_cut_hpf = 1.6e6; + +[zb1, pb1, kb1] = butter(1, 2*pi*f_cut_hpf, 'high', 's'); +[zd1, pd1, kd1] = bilinear(zb1, pb1, kb1, 2*fs); +[hpf1_d_num, hpf1_d_den] = zp2tf(zd1, pd1, kd1); + +temp_mixed_sig_hpf = zeros(size(mixed_sig)); +mixed_sig_hpf = zeros(size(mixed_sig)); +temp_mixed_noise_hpf = zeros(size(mixed_sig)); +mixed_noise_hpf = zeros(size(mixed_sig)); +for rxidx = 1 : Nr + if opt_filter == 1 + temp_mixed_sig_hpf(:,:,rxidx) = filter(hpf1_d_num, hpf1_d_den, mixed_sig(:,:,rxidx)); + mixed_sig_hpf(:,:,rxidx) = filter(hpf1_d_num, hpf1_d_den, temp_mixed_sig_hpf(:,:,rxidx)); + + temp_mixed_noise_hpf(:,:,rxidx) = filter(hpf1_d_num, hpf1_d_den, mixed_noise(:,:,rxidx)); + mixed_noise_hpf(:,:,rxidx) = filter(hpf1_d_num, hpf1_d_den, temp_mixed_noise_hpf(:,:,rxidx)); + else + mixed_sig_hpf(:,:,rxidx) = mixed_sig(:,:,rxidx); + + mixed_noise_hpf(:,:,rxidx) = mixed_noise(:,:,rxidx); + end +end +[h_hpf1, w] = freqz(hpf1_d_num,hpf1_d_den, rngFFTlen); +Loss_hpf1_dB = pow2db(sum(abs(h_hpf1).^2)/rngFFTlen); +Loss_hpf2_dB = pow2db(sum(abs(h_hpf1).^2)/rngFFTlen); + + + +% 2-3. Low Pass Filter +% 1) Third-order, -6 dB, (12.5 MHz ~ 25 MHz) : NXP (tef82xx) +% 2) Third-order, -6 dB, > 40 MHz(Wide-bandwidth mode) : NXP (tef82xx) +f_cut_lpf = 12.5e6; % (20e6 * 416/512) + +[zl1, pl1, kl1] = butter(1, 2*pi*f_cut_lpf, 'low', 's'); +[zdl1, pdl1, kdl1] = bilinear(zl1, pl1, kl1, 2*fs); +[lpf1_d_num, lpf1_d_den] = zp2tf(zdl1, pdl1, kdl1); + +[zl2, pl2, kl2] = butter(2, 2*pi*f_cut_lpf, 'low', 's'); +[zdl2, pdl2, kdl2] = bilinear(zl2, pl2, kl2, 2*fs); +[lpf2_d_num, lpf2_d_den] = zp2tf(zdl2, pdl2, kdl2); + +temp_mixed_sig_lpf = zeros(size(mixed_sig)); +mixed_sig_lpf = zeros(size(mixed_sig)); +temp_mixed_noise_lpf = zeros(size(mixed_sig)); +mixed_noise_lpf = zeros(size(mixed_sig)); +for rxidx = 1 : Nr + if opt_filter == 1 + temp_mixed_sig_lpf(:,:,rxidx) = filter(lpf1_d_num, lpf1_d_den, mixed_sig_hpf(:,:,rxidx)); + mixed_sig_lpf(:,:,rxidx) = filter(lpf2_d_num, lpf2_d_den, temp_mixed_sig_lpf(:,:,rxidx)); + + temp_mixed_noise_lpf(:,:,rxidx) = filter(lpf1_d_num, lpf1_d_den, mixed_noise_hpf(:,:,rxidx)); + mixed_noise_lpf(:,:,rxidx) = filter(lpf2_d_num, lpf2_d_den, temp_mixed_noise_lpf(:,:,rxidx)); + else + mixed_sig_lpf(:,:,rxidx) = mixed_sig_hpf(:,:,rxidx); + + mixed_noise_lpf(:,:,rxidx) = mixed_noise_hpf(:,:,rxidx); + end +end +[h_lpf1, w] = freqz(lpf1_d_num, lpf1_d_den, rngFFTlen); +Loss_lpf1_dB = pow2db(sum(abs(h_lpf1).^2)/rngFFTlen); +[h_lpf2, w] = freqz(lpf2_d_num, lpf2_d_den, rngFFTlen); +Loss_lpf2_dB = pow2db(sum(abs(h_hpf1).^2)/rngFFTlen); + +% figure() +% subplot(211); plot(Tx_wave.timeline, real(mixed_sig_lpf(:,10))); +% xlabel('Time (s)'); ylabel('Amplitude (v)'); +% title('dechirped FMCW signal'); axis tight; +% subplot(212); spectrogram(mixed_sig_lpf(:,10), 320, 160, 320, fs,'yaxis'); +% title('dechirped FMCW signal spectrogram'); + +% 2-4. Receiver gain +Rxgain_dB = 45; +mixed_sig_lpf_rxgain = mixed_sig_lpf(1:2:end,:,:) * db2pow((Rxgain_dB - (Loss_hpf1_dB + Loss_hpf1_dB + Loss_lpf1_dB + Loss_lpf2_dB))/2); +mixed_noise_lpf_rxgain = mixed_noise_lpf(1:2:end,:,:) * db2pow((Rxgain_dB - (Loss_hpf1_dB + Loss_hpf1_dB + Loss_lpf1_dB + Loss_lpf2_dB))/2); + +% 2-5 Quantization Noise +Vpp_ADC = 1.2; +ADC_bit = 12; +ADC_impd = 50; +Vq_rms = Vpp_ADC/2^ADC_bit / sqrt(12); +qt_pow = Vq_rms^2 / ADC_impd; +qt_noise = sqrt(qt_pow) * randn(size(mixed_sig_lpf_rxgain)); + +sig_adc_out = mixed_sig_lpf_rxgain; +noise_adc_out = mixed_noise_lpf_rxgain; + +%% 03. FFT +% 3-1 Range Windowing +if opt_window == 1 + winRng = repmat(hann(rngFFTlen), [1,NumChirps, Nr]); +else + winRng = repmat(ones(NumSamples, 1), [1,NumChirps, Nr]); +end + +scalwinRng = sum(winRng(:,1,1)) / length(winRng(:,1,1)); +winR_mixed_sig_lpf = sig_adc_out .* winRng; +winR_mixed_noise_lpf = noise_adc_out .* winRng; + +% 3-2. Range FFT +sig_rng_fft = fft(winR_mixed_sig_lpf, rngFFTlen, 1); +noise_rng_fft = fft(winR_mixed_noise_lpf, rngFFTlen, 1); + +% 3-3 Doppler Windowing +if opt_window == 1 + winDop = repmat(hann(velFFTlen).', [rngFFTlen, 1, Nr]); +else + winDop = repmat(ones(1, NumChirps), [rngFFTlen, 1, Nr]); +end + +scalwinDop = sum(winDop(1,:,1)) / length(winDop(1,:,1)); +winD_sig_rng_fft = sig_rng_fft .* winDop; +winD_noise_rng_fft = noise_rng_fft .* winDop; + +% 3-4. Doppler FFT +sig_rng_dop_fft2_wo_n = fft(winD_sig_rng_fft, velFFTlen, 2); +noise_rng_dop_fft2 = fft(winD_noise_rng_fft, velFFTlen, 2); + +%noise_att =[0.2120 0.2143 0.2206 0.2276 0.2342 0.2424 0.2468 ]; +noise_att = ones(1, maxrngFFTidx); +noise_rng_dop_fft2(1:maxrngFFTidx, :, :) = noise_rng_dop_fft2(1:maxrngFFTidx, : ,:) .* repmat(noise_att', 1, 256, 4); + +sig_rng_dop_fft2 = sig_rng_dop_fft2_wo_n + noise_rng_dop_fft2; + + + + +%% Signal Processing +rng_grid = beat2rng(gen_freqgrid(rngFFTlen, fs, 0), Tx_wave.f_slope, c0); +spd_grid = 0.5*dop2spd(gen_freqgrid(velFFTlen, 1/(Tx_wave.T_chirp), 1), lambda_c); +spd_grid_2 = fftshift(spd_grid); + +PS_sig_single_ch = abs(fftshift(sig_rng_dop_fft2 ,2)).^2/NumChirps/NumSamples/rngFFTlen/velFFTlen; +PS_noise_single_ch = abs(fftshift(noise_rng_dop_fft2 ,2)).^2/NumChirps/NumSamples/rngFFTlen/velFFTlen; + +Winloss_dB = 0; +RBW_noise = fs/velFFTlen/rngFFTlen; +%noise_floor_dBm = pow2db(kb * T0) + NF_dB + pow2db(RBW_noise) + Rxgain_dB + Winloss_dB + pow2db(sum(winDop(1,:,1).^2)/velFFTlen) + pow2db(sum(winRng(:,1,1).^2)/rngFFTlen) + 30; +%noise_floor_dBm = pow2db(kb * T0) + NF_dB + pow2db(RBW_noise) + Rxgain_dB + Winloss_dB + pow2db(1/4) - 1.76*2 + 30; +noise_floor_dBm = pow2db(kb * T0) + NF_dB + pow2db(RBW_noise) + Rxgain_dB + pow2db(sum(winDop(1,:,1).^2)/velFFTlen) + pow2db(sum(winRng(:,1,1).^2)/rngFFTlen) + 30; + +% Mean noise power in single channel +np_true_single_ch = sum(PS_noise_single_ch,2) / NumChirps; + +np_esti_single_ch = zeros(maxrngFFTidx, Nr); +for ch_idx = 1 : Nr + for rng_idx = 1 : maxrngFFTidx + [find_noise_index, l, u, np_esti_single_ch(rng_idx, ch_idx)] = isoutlier(PS_sig_single_ch(rng_idx, :, ch_idx), "percentiles", [10 80]); + end +end + + +if flag_plot_RVM == 1 + figure + for ch_idx = 1 : Nr + subplot(sqrt(Nr), sqrt(Nr), ch_idx); + mesh(spd_grid, rng_grid(1:maxrngFFTidx), pow2db(squeeze(PS_sig_single_ch(1:maxrngFFTidx, :, ch_idx))) + 30); + xlabel('Rng [m]'); + ylabel('Vel [m/s]'); + zlabel('Pow [dBm]'); + hold on + mesh(spd_grid, rng_grid(1:maxrngFFTidx), noise_floor_dBm*ones(maxrngFFTidx, velFFTlen), 'EdgeColor', 'r'); + mesh(spd_grid, rng_grid(1:maxrngFFTidx), pow2db(np_esti_single_ch(:, ch_idx).*ones(1,velFFTlen)) + 30, 'EdgeColor', 'g'); + legend('Received Data', 'Expected Noise Floor', 'Esti. Noise Floor', 'Location', 'best'); + hold off + title(['Ch : ', num2str(ch_idx)]); + end +end + +% figure +% imagesc(spd_grid, rng_grid(1:416), (pow2db(abs(sig_fft2(1:416,:,1))))) +% xlabel('Speed (m/s)'); ylabel('Range (m)'); title('Range Velocity Map'); +% axis([-vel_max vel_max 0 rng_max]) +% colorbar + +% 1. Generate RV Quarter Matrix +% 1-1. NCI Rx +sig_Rx_NCI = sum(abs(sig_rng_dop_fft2).^2, 3) / Nr; +noise_Rx_NCI = sum(abs(noise_rng_dop_fft2).^2, 3) / Nr; + +PS_sig_Rx_NCI = sig_Rx_NCI/NumChirps/NumSamples/rngFFTlen/velFFTlen; +PS_noise_Rx_NCI = noise_Rx_NCI/NumChirps/NumSamples/rngFFTlen/velFFTlen; + +if flag_plot_RVM == 1 + figure + mesh(spd_grid, rng_grid(1:maxrngFFTidx), pow2db(fftshift(PS_sig_Rx_NCI(1:maxrngFFTidx,:), 2)) + 30); + axis([-vel_max vel_max 0 rng_max]); + xlabel('Speed [m/s]'); + ylabel('Range [m]'); + zlabel('Power [dBm]'); + title('RVM - NCI w.r.t Rx'); +end + +% NCI Tx +sig_TxRx_NCI = (1/4) * (sig_Rx_NCI(:,1:velFFTlen/4) + sig_Rx_NCI(:,velFFTlen/4 + 1 : 2*velFFTlen/4) + sig_Rx_NCI(:,2*velFFTlen/4 + 1 : 3*velFFTlen/4 ) + sig_Rx_NCI(:,3*velFFTlen/4+1 : end)); +noise_TxRx_NCI = (1/4) * (noise_Rx_NCI(:,1:velFFTlen/4) + noise_Rx_NCI(:,velFFTlen/4 + 1 : 2*velFFTlen/4) + noise_Rx_NCI(:,2*velFFTlen/4 + 1 : 3*velFFTlen/4 ) + noise_Rx_NCI(:,3*velFFTlen/4+1 : end)); + +PS_sig_TxRx_NCI = sig_TxRx_NCI/NumChirps/NumSamples/rngFFTlen/velFFTlen; +PS_noise_TxRx_NCI = noise_TxRx_NCI/NumChirps/NumSamples/rngFFTlen/velFFTlen; + +if flag_plot_RVM == 1 + figure + mesh(pow2db(PS_sig_TxRx_NCI(1:maxrngFFTidx,:)) + 30); + xlabel('Dop idx'); + ylabel('Rng idx'); + zlabel('Power [dBm]'); + title('RVM - NCI w.r.t Tx and Rx'); +end + +% 2. Noise Estimation +np_esti_TxRx_NCI = zeros(maxrngFFTidx, 1); +for rng_idx = 1 : maxrngFFTidx + [cal_noise_index, l, u, np_esti_TxRx_NCI(rng_idx)] = isoutlier(PS_sig_TxRx_NCI(rng_idx, :), "percentiles", [10 60]); +end + +np_true_TxRx_NCI = sum(PS_noise_TxRx_NCI(1:maxrngFFTidx,:), 2) / (NumChirps/(length(DDMA_idx)+1)); + +if flag_plot_noise_esti == 1 + figure + plot(pow2db(np_true_single_ch(1:maxrngFFTidx)) + 30); + hold on + plot(pow2db(np_esti_single_ch) + 30); + plot(pow2db(np_true_TxRx_NCI(1:maxrngFFTidx)) + 30); + plot(pow2db(np_esti_TxRx_NCI) + 30); + legend('NP ture in single ch.', 'NP esti in single ch.', 'NP true at TxRxNCI', 'NP esti at TxRxNCI'); + zlabel('Power [dBm]'); +end + +% 3. Detection +% 3-1. Setting Threshold +snr_th = db2pow(13); + +% 3-2. Find det index +[det_index_peak, b] = findpeak2D(PS_sig_TxRx_NCI(1:maxrngFFTidx,:), []); + +det_index_th = []; +for rngidx = 1 : maxrngFFTidx + det_index_th = [det_index_th ; PS_sig_TxRx_NCI(rngidx, :) > np_esti_TxRx_NCI(rngidx) * snr_th]; +end + +det_index = det_index_th .* det_index_peak; + +[rngidx_set, dopidx_set] = find(det_index>0); + +for det_idx = 1 : length(rngidx_set) + temp_DET_set{det_idx}.r_m = rng_grid(rngidx_set(det_idx)); + temp_DET_set{det_idx}.snr_db = pow2db((PS_sig_TxRx_NCI(rngidx_set(det_idx), dopidx_set(det_idx))) / np_esti_TxRx_NCI(rngidx_set(det_idx))); +end + +% 4. DOA estimation +% 4-1. Generate Snapshot and Resolving Doppler ambiguity by DDMA +adddopidx = [64 128 192 0]; +snapshot_idx_set = [4 1 2 3; 3 4 1 2; 2 3 4 1; 1 2 3 4]; +temp_snapshot_data = zeros((Nt+1)*Nr,length(rngidx_set)); +for idx = 1 : length(rngidx_set) + temp_snapshot = zeros((Nt+1)*Nr,1); + for txidx = 1 : Nt+1 + temp_snapshot(Nr*(txidx-1)+1 : Nr*txidx) = squeeze(sig_rng_dop_fft2_wo_n(rngidx_set(idx), dopidx_set(idx)+64*(txidx-1), :)); + end + + % Find virtual channel by selecting min power channel + rxchpow = sum(abs(reshape(temp_snapshot, Nt+1, Nr)).^2); + [minval, minloc] = min(rxchpow); + + dopidx_set(idx) = dopidx_set(idx) + adddopidx(minloc); + + for txidx = 1 : Nt+1 + temp_snapshot_data(Nr*(txidx-1)+1 : Nr*txidx, idx) = temp_snapshot(4*snapshot_idx_set(minloc, txidx)-3 : 4*snapshot_idx_set(minloc, txidx)); + end +end +snapshot_data = temp_snapshot_data(1:Nt*Nr, :); + +% Add resolved vel info to DET_set structure +for det_idx = 1 : length(rngidx_set) + temp_DET_set{det_idx}.v_amb_mps = spd_grid_2(dopidx_set(det_idx)); +end + +% 4-2. DOA estimation +az_array_pos = mimoarray(2, :) / u_azi; + +elv_Unit = u_elv / lambda_c; + +test_ang = -90 : 0.1 : 90; +svmat = exp(-1i * 2 * pi / lambda_c * az_array_pos.' * u_azi .* sind(test_ang)); + +det_jdx = 1; +for det_idx = 1 : size(snapshot_data, 2) + % DOA estimation + % Elv esti. + elv_phase_diff = conj(snapshot_data(4, det_idx)) * snapshot_data(9, det_idx); + temp_esti_elv_deg = asind(angle(elv_phase_diff)/2/pi/elv_Unit); + + % Azi esti. + azi_spectrum = abs(svmat' * snapshot_data(:, det_idx)).^2; + [pks, esti_azi_deg] = findpeaks(azi_spectrum/max(azi_spectrum), test_ang, 'MinPeakHeight', 0.5); + + esti_elv_deg = temp_esti_elv_deg * ones(1, length(esti_azi_deg)); + + + % Construct DET set structure + for angidx = 1 : length(esti_azi_deg) + DET_set{det_jdx} = temp_DET_set{det_idx}; + DET_set{det_jdx}.az_deg = esti_azi_deg(angidx); + DET_set{det_jdx}.el_deg = esti_elv_deg(angidx); + DET_set{det_jdx}.x_m = DET_set{det_jdx}.r_m * cosd(DET_set{det_jdx}.az_deg); + DET_set{det_jdx}.y_m = DET_set{det_jdx}.r_m * sind(DET_set{det_jdx}.az_deg); + DET_set{det_jdx}.z_m = DET_set{det_jdx}.r_m * sind(DET_set{det_jdx}.el_deg); + det_jdx = det_jdx + 1; + end +end + +%% Results +if flag_plot_birdview == 1 + figure; + x_m_set = cellfun(@(x) x.x_m, DET_set); + y_m_set = cellfun(@(x) x.y_m, DET_set); + z_m_set = cellfun(@(x) x.z_m, DET_set); + r_m_set = cellfun(@(x) x.r_m, DET_set); + v_amb_mps_set = cellfun(@(x) x.v_amb_mps, DET_set); + az_deg_set = cellfun(@(x) x.az_deg, DET_set); + el_deg_set = cellfun(@(x) x.el_deg, DET_set); + snr_db_set = cellfun(@(x) x.snr_db, DET_set); + plotdet = plot3(y_m_set, x_m_set, z_m_set, 'bo'); + fn_add_data_tip(plotdet, 'Rng [m]:', r_m_set, 1); + fn_add_data_tip(plotdet, 'aVel [m/s]:', v_amb_mps_set, 2); + fn_add_data_tip(plotdet, 'Azi [deg]:', az_deg_set, 3); + fn_add_data_tip(plotdet, 'Elv [deg]:', el_deg_set, 4); + fn_add_data_tip(plotdet, 'SNR [dB]:', snr_db_set, 5); + fn_add_data_tip(plotdet, 'X [m]:', x_m_set, 6); + fn_add_data_tip(plotdet, 'Y [m]:', y_m_set, 7); + fn_add_data_tip(plotdet, 'Z [m]:', z_m_set, 8); + view([0, 90]); + grid on + xlabel('Y [m]') + ylabel('X [m]') + zlabel('Z [m]') + xlim([-rng_max, rng_max]) + ylim([0 rng_max]) + title('BirdView') +end + + +