function [PDP, rmsDelay, psd_max] = rxChannelSounding(rx,cfg,numSeq,K)
chanBW = cfg.ChannelBandwidth;
T_sample = 1/str2double(chanBW(4:end))/1e6; % Duration of each sample/pulse
% K = 255; % Length of a pseudonoise sequence
T_seq = K*T_sample; % Duration of a pseudonoise sequence

soundingSeqMod_1 = zadoffChuSeq(29,255);
soundingSeqMod_2 = zadoffChuSeq(47,255);
%% 计算PDP
% Concatenate multiple sequences to sound the channel periodically.
% Preprocessing of the Received Data
MISO_channel_index = 1 ;%/ 2
if MISO_channel_index == 1
    soundingSeqMod = soundingSeqMod_1;
else
    soundingSeqMod = soundingSeqMod_2;
end

% Create a matched filter of the modulated sounding sequence.
h_mf = conj(flip(soundingSeqMod));

% Pass the received signal through the matched filter
output_mf = conv(h_mf,rx)/K;
output_mf = output_mf(K:end);

% Store the real part and imaginary part of the received data respectively.
data_I = real(output_mf);              % I route for recorded data
data_Q = imag(output_mf);              % Q route for recorded data
% figure;
% Power Delay Profiles (PDPs)
PDP = mean(reshape(data_I.^2+data_Q.^2,K,numSeq),2);
x = (0:K-1)*T_sample*1e6;
% stem(x,PDP);
% axis([0 T_seq*1e5 0 1.2*max(PDP)]);
% xlabel('delay (us)');
% ylabel('power (W)');
PDP2d = reshape(data_I.^2+data_Q.^2,K,numSeq);
% pause(0.5);
% close;
% [x,y] = meshgrid(1:numSeq,(0:K-1)*T_sample*1e6);
% mesh(x,y,PDP)
% axis([1 numSeq 0 T_seq*1e5 0 1.5]);
% xlabel('Index of PDP')
% ylabel('delay (us)')
% zlabel('power (W)')
% title('Power Delay Profiles (PDP)');

%% 计算时延扩展
% Compute the temporal parameters.
delay_average_vector = zeros(numSeq,1);
delay_spread_vector = zeros(numSeq,1);

for i = 1:numSeq
    [delay_average_vector(i), delay_spread_vector(i)] = findTemporalParameters(PDP2d(:,i), T_sample, K);
end
range_delay = 0:T_sample:T_seq;
% pdf_delay_average = findPDF(range_delay,delay_average_vector);
pdf_delay_spread = findPDF(range_delay,delay_spread_vector);
rmsDelay = range_delay(pdf_delay_spread==max(pdf_delay_spread))*1e6;
% figure;
% subplot(1,2,1)
% plot(range_delay(1:length(range_delay)-1)*1e6,pdf_delay_average,'LineWidth',2);
% axis([0 1.25*max(delay_average_vector)*1e6 0 1.25*max(pdf_delay_average)])
% xlabel('Delay (us)')
% ylabel('Average Delay')
% grid on
% subplot(1,2,2)
% plot(range_delay(1:length(range_delay)-1)*1e6,pdf_delay_spread,'LineWidth',2);
% xlabel('Delay (us)')
% ylabel('Delay Spread')
% axis([0 1.25*max(delay_spread_vector)*1e6 0 1.25*max(pdf_delay_spread)])
% grid on

%% 计算多普勒频移
% Truncate the measurement data into individual impulse response for further analysis
dataBlk = reshape(data_I+sqrt(-1)*data_Q,K,numSeq);

% Compute the Doppler power spectral density through FFT of the ACF.
% What if the delay of particular path varies with time?
psd = zeros(K,2*numSeq-1);                  % Doppler PSD
for i = 1:K
    temp = xcorr(dataBlk(i,:));
    psd(i,:) = fftshift(abs(fft(temp)/(2*numSeq-1)));
end
f = ((0:2*numSeq-2)/(2*numSeq-1)-0.5)/T_seq;
psd_max = f(max(psd)==max(psd,[],'all'));
% delay = [0 0.2 0.4 0.6]*1e-6;
% figure;
% [x,y] = meshgrid(((0:2*numSeq-2)/(2*numSeq-1)-0.5)/T_seq,(0:K-1)*T_sample*1e6);
% mesh(x,y,psd)
% axis([-1000 1000 0 1.5*max(delay)*1e6 0 1.5*max(max(psd))])
% xlabel('Doppler Shift (Hz)')
% ylabel('Delay (us)')
% zlabel('Doppler Power Spectral Density')
end