Initial Commit.

2025-07-29 21:45:44 +09:00
parent a69350b86e
commit 90231d49ae
38 changed files with 2051 additions and 0 deletions
@@ -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
+
+
+
+
+
@@ -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)
@@ -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
+
@@ -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
@@ -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
+
@@ -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
@@ -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
@@ -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
+
@@ -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
+
@@ -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
+
@@ -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
@@ -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
+
@@ -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
+
@@ -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
+
@@ -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
@@ -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
@@ -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
+
@@ -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
+
@@ -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
+
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
+
@@ -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
+
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
+
@@ -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
+