clc; clear all; close all;
load('bellhop_dataset.mat'); % 加载扫频数据集，dataset是结构体数组 扫频步长frequencies = 3500:100:6500;
Fc = 5000; % 中心频率
bw = 2000; % 带宽
% 定义频率区间（20个，步长100 Hz，覆盖Fc - bw/2 到 Fc + bw/2）
interval_size = 100; % 区间大小 (Hz)
num_intervals = bw / interval_size;
interval_centers = Fc - bw/2 + interval_size/2 : interval_size : Fc + bw/2 - interval_size/2; % 中心频率
Rs = 500; % 根据香农定理 此时符号带宽为500波特
V_rate = 1000; % 通讯速率
Fs = 16000; % 采样率
sl_db = 157; % 发射机声源级
windspeed = 5; % 海面风浪等级
shipping_factor = 0.5;% 海面船只密度 用于评估船噪声
c_water = 1500; % 水声速 m/s
% 发射端和接收端速度向量 空间向量
V_tx = [0, 0, 0]; % 发射端速度向量 [vx, vy, vz] m/s（正：向接收端接近）
V_rx = [0, 0, 0]; % 接收端速度向量 [vx, vy, vz] m/s（正：向发射端接近）
V_current = [0, 0, 0]; % 洋流速度向量 m/s
% OFDM 参数
para = 1024; % 子载波数
fftlen = 1024; % FFT长度
gilen = 1024; % 保护间隔长度
nd = 80; % 每循环OFDM符号数
ml = 2; % QPSK调制，2比特/符号
fftlen2 = fftlen + gilen; % 总符号长度
% 添加导频参数
pilot_spacing = 4; % 导频间隔（一般4~16）
num_pilots = ceil(para / pilot_spacing); % 导频数
pilot_indices = 1:pilot_spacing:para; % 导频位置（从1开始，每间隔一个）
pilot_value = 1 / sqrt(2); % 导频值（功率归一化，与QPSK一致）
% 添加null子载波参数（用于噪声估计）
null_spacing = 16; % 避免与pilot重叠
num_nulls = ceil(para / null_spacing); %空载波数
null_indices = pilot_spacing/2 : null_spacing : para; %空载波索引值 从中间起始，避免重叠（e.g., 4,20,36,...）
% 噪声估计选项（为脉冲噪声添加median）
use_median_for_noise = false; % 如果您的噪声有脉冲，设为true；否则false用mean
% 验证导频与空载波无重叠
if ~isempty(intersect(pilot_indices, null_indices))
error('导频索引与零值索引存在重叠！请调整起始位置');
end
% 更新数据子载波
data_subcarriers = setdiff(1:para, union(pilot_indices, null_indices));
num_data_sc = length(data_subcarriers);
% 计算OFDM符号时长和sps（基于带宽bw）
subcarrier_spacing = bw / para; % 子载波间隔 (Hz)
T_eff = fftlen / bw; % 有效符号时长 (s)
T_gi = gilen / bw; % GI时长 (s)
T_symbol = T_eff + T_gi; % 总符号时长 (s)
sps = round(Fs * T_symbol); % 每个OFDM符号采样数
% 深度估计（全局，从第一个取）
position_s = Belldataset(1).arrivals.Pos.s.z;
position_r = Belldataset(1).arrivals.Pos.r.z;
position_rx = Belldataset(1).arrivals.Pos.r.r;
depth = ( position_r + position_s )/2; % 粗略估计损耗位置 单位m
num_paths = Belldataset(1).arrivals.Arr.Narr;
% 预分配预处理结构体数组（每个区间一个）
preprocessed_data = struct('sorted_A', cell(1, num_intervals), ...
'norm_delay', cell(1, num_intervals), ...
'sorted_SrcAzimuthAngle', cell(1, num_intervals), ...
'sorted_RcvrAzimuthAngle', cell(1, num_intervals), ...
'fd_paths', cell(1, num_intervals),...
'num_paths', zeros(1, num_intervals));
for int_idx = 1:num_intervals
% 检索最近的
[~, data_idx] = min(abs([Belldataset.frequency] - interval_centers(int_idx)));
% 检查越界
if data_idx > length(Belldataset)
error('data_idx %d 超过 Belldataset 长度 %d，请调整偏移', data_idx, length(Belldataset));
end
% 提取数据
Arr = Belldataset(data_idx).arrivals.Arr;
Pos = Belldataset(data_idx).arrivals.Pos; % 备选
A_complex = Arr.A; % 复振幅
delay = Arr.delay; % 时延
SrcAzimuthAngle = Arr.SrcDeclAngle; % 出射角
RcvrAzimuthAngle = Arr.RcvrDeclAngle; % 入射角
% 升序排序
[~, sort_idx] = sort(delay);
sorted_delay = delay(sort_idx); % 排序后的时间
sorted_A = A_complex(sort_idx); % 排序后的幅值
sorted_SrcAzimuthAngle = SrcAzimuthAngle(sort_idx); % 排序后的出射角
sorted_RcvrAzimuthAngle = RcvrAzimuthAngle(sort_idx); % 排序后的入射角
norm_delay = sorted_delay - sorted_delay(1); % 归一化延时（假设至少一条路径）
% 计算 per区间 多普勒（用区间中心频率）
num_paths = length(norm_delay);
center_freq = interval_centers(int_idx); % 用区间中心频替换 Fc
fd_paths = zeros(num_paths, 1);
for k = 1:num_paths
theta_s_deg = sorted_SrcAzimuthAngle(k);
theta_r_deg = sorted_RcvrAzimuthAngle(k); % 角度获取
theta_s = deg2rad(theta_s_deg);
theta_r = deg2rad(theta_r_deg); % 角度转换
theta_avg = (theta_s + theta_r) / 2; % 粗略估计
u_k = [cos(theta_avg), 0, sin(theta_avg)]; % 方向向量
v_s_radial = dot(V_tx, u_k);
v_r_radial = dot(V_rx, u_k);
v_c_radial = dot(V_current, u_k); % 点积
fd_paths(k) = (v_s_radial + v_r_radial + v_c_radial) / c_water * center_freq; % 总多普勒频移，用 center_freq
end
% 存储
preprocessed_data(int_idx).sorted_A = sorted_A;
preprocessed_data(int_idx).norm_delay = norm_delay;
preprocessed_data(int_idx).sorted_SrcAzimuthAngle = sorted_SrcAzimuthAngle;
preprocessed_data(int_idx).sorted_RcvrAzimuthAngle = sorted_RcvrAzimuthAngle;
preprocessed_data(int_idx).fd_paths = fd_paths; % 新增存储
preprocessed_data(int_idx).num_paths = num_paths;
% 进度（可选，新增显示 fd_paths 范围验证）
% fprintf('预处理完成区间 %d (中心频 %.0f Hz, dataset索引 %d, 频 %.0f Hz), 路径数: %d, 多普勒范围: [%.3f, %.3f] Hz\n', ...
% int_idx, center_freq, data_idx, 3500 + 100*(data_idx-1), num_paths, min(fd_paths), max(fd_paths));
end
% 在OFDM参数定义后，预计算CIR长度L（全局，稀疏估计用）
max_delay_all = 0;
for int_idx = 1:num_intervals
max_delay_all = max(max_delay_all, max(preprocessed_data(int_idx).norm_delay));
end
L = round(max_delay_all * bw * 1.2); % CIR长度，裕度1.2倍（基于带宽bw，采样到基带率）
% % 在预处理循环前添加（覆盖所有区间）
% for int_idx = 1:num_intervals
% preprocessed_data(int_idx).sorted_A = 0.000001; % 单径，幅度1（无损耗）
% preprocessed_data(int_idx).norm_delay = 0;
% preprocessed_data(int_idx).fd_paths = 0;
% preprocessed_data(int_idx).num_paths = 1;
% % 其他字段保持
% end
%%
% === 3. 数据生成 ===
% 说明：生成随机比特，长度适配OFDM：总比特数 = 子载波数 × OFDM符号数 × 每符号比特数。
% 原因：OFDM将数据分配到多个子载波并行传输，需计算总比特。
N_bits = num_data_sc * nd * ml;
source = randsrc(1, N_bits, [1, 0]); % 随机二进制序列
% === 4. 发射机：OFDM调制 ===
% 说明：实现串并转换、QPSK调制、IFFT、添加保护间隔
% 原因：OFDM发射机将比特映射到频域子载波，转换为时域信号，添加保护间隔对抗多径。
paradata = reshape(source, num_data_sc, nd*ml); % 串并转换：XX行，20列
% QPSK调制
[ich, qch] = qpskmod(paradata, num_data_sc, nd, ml);
kmod = 1/sqrt(2); % 归一化因子
ich1 = ich * kmod;
qch1 = qch * kmod;
x_data = ich1 + qch1 * 1j; % 频域复信号 数据部分：num_data_sc × nd
% 插入导频和空载波：创建全x矩阵
x = zeros(para, nd); % 全频域信号
x(data_subcarriers, :) = x_data; % 填入数据
x(pilot_indices, :) = pilot_value; % 填入导频（所有符号相同）
x(null_indices, :) = 0; % null设为零
% IFFT
y = ifft(x, fftlen) * sqrt(fftlen); % 修改：归一ifft
ich2 = real(y); %实数时域信号
qch2 = imag(y); %虚部时域信号
% 添加保护间隔
[ich3, qch3] = giins(ich2, qch2, fftlen, gilen, nd); % 输出：80×10
% 转换为单载波时域信号并上变频
tx_sig = zeros(1, nd * sps);
t_symbol = (0:sps-1)/Fs;
t_base = (0:fftlen2-1)/bw;
for i = 1:nd
start_idx = (i-1)*sps + 1;
end_idx = i*sps;
symbol = ich3(:, i) + 1j * qch3(:, i); % 基带符号，长80
% 计算上采样因子
up_factor = Fs / bw;
% resample上采样基带复信号（symbol是复数）
symbol_interp = resample(symbol, up_factor, 1);
% 确保长度精确（resample可能因滤波略微调整，trim或pad如果需要）
if length(symbol_interp) > sps
symbol_interp = symbol_interp(1:sps);
elseif length(symbol_interp) < sps
symbol_interp = [symbol_interp; zeros(sps - length(symbol_interp), 1)];
end
% 上变频并取实部（t_symbol不变）
symbol_rf = real(symbol_interp .' .* exp(1j * 2 * pi * Fc * t_symbol)); % 注意转置为行向量
tx_sig(start_idx:end_idx) = symbol_rf;
end
% === 5. 信道：多径和多普勒（宽带分段版本） ===
% 先计算全局 max_delay_all（所有区间最大，用于缓冲）
max_delay_all = 0;
for int_idx = 1:num_intervals
max_delay_all = max(max_delay_all, max(preprocessed_data(int_idx).norm_delay));
end
max_delay_samples = round(max_delay_all * Fs);
padding_len = max_delay_samples; % 定义 padding 长度（保守取 max_delay_samples）
% 填充 tx_sig 后加零
tx_sig_padded = [tx_sig, zeros(1, padding_len)]; %填充零
padded_len = length(tx_sig_padded); % 填充后长度
% 将 tx_sig_padded 转换为频域（用 padded_len）
N_fft = padded_len; % 用填充后长度
tx_freq = fft(tx_sig_padded, N_fft); % 用填充信号 FFT
tx_freq_shifted = fftshift(tx_freq); % 保持
f_axis = (0:N_fft-1) * (Fs / N_fft) - Fs/2; % 基于新 N_fft
% 为每个区间构建 H(f) 并应用
tx_freq_channel_shifted = zeros(size(tx_freq_shifted));
for int_idx = 1:num_intervals
center_freq = interval_centers(int_idx);
f_low = center_freq - interval_size/2;
f_high = center_freq + interval_size/2;
% 从预处理调用
sorted_A = preprocessed_data(int_idx).sorted_A;
norm_delay = preprocessed_data(int_idx).norm_delay;
fd_paths = preprocessed_data(int_idx).fd_paths;
num_paths = preprocessed_data(int_idx).num_paths;
% 构建区间 H(f)
bin_mask_pos = f_axis >= f_low & f_axis < f_high; % 正频
bin_mask_neg = f_axis <= -f_low & f_axis > -f_high; % 负频
f_bins_pos = f_axis(bin_mask_pos);
f_bins_neg = f_axis(bin_mask_neg);
% 正频 H_f_pos
H_f_pos = zeros(size(f_bins_pos)) + 1j * zeros(size(f_bins_pos));
for bin_idx = 1:length(f_bins_pos)
f = abs(f_bins_pos(bin_idx)); % 正频abs
for k = 1:num_paths
A_k = sorted_A(k) * 10^(sl_db / 20);
tau_k = norm_delay(k);
fd_k = fd_paths(k);
H_f_pos(bin_idx) = H_f_pos(bin_idx) + A_k * exp(-1j * 2 * pi * (f + fd_k) * tau_k);
end
end
% 负频 H_f_neg
H_f_neg = zeros(size(f_bins_neg)) + 1j * zeros(size(f_bins_neg));
for bin_idx = 1:length(f_bins_neg)
f = f_bins_neg(bin_idx); % 用负值
for k = 1:num_paths
A_k = sorted_A(k) * 10^(sl_db / 20);
tau_k = norm_delay(k);
fd_k = fd_paths(k);
H_f_neg(bin_idx) = H_f_neg(bin_idx) + A_k * exp(-1j * 2 * pi * (f + fd_k) * tau_k);
end
end
tx_freq_channel_shifted(bin_mask_pos) = tx_freq_shifted(bin_mask_pos) .* H_f_pos;
tx_freq_channel_shifted(bin_mask_neg) = tx_freq_shifted(bin_mask_neg) .* H_f_neg;
end
% IFFT 回时域，取实部
tx_freq_channel = ifftshift(tx_freq_channel_shifted); % 保持
tx_sig_channel_full = real(ifft(tx_freq_channel, N_fft)); % 修改：用 N_fft IFFT，得填充长输出
% 截取：线性卷积输出长 = 原长 + padding_len，但我们只需前 length(tx_sig) + max_delay_samples - 1
tx_sig_channel = tx_sig_channel_full(1:length(tx_sig)); % 修改：截取原长（保守）
% === 6. 噪声处理 ===
% per区间 平均噪声 PSD 和 maxamp_db/ir_vec（捕捉频率变异）
npsd_db_avg = 0;
maxamp_db = -Inf; % 初始化为最小
ir_vec_cell = cell(1, num_intervals); % 用 cell 收集列向量
max_ir_len = 0; % 找最大长度
for int_idx = 1:num_intervals
[npsd_db_int] = improved_ambientnoise_psd_v2(windspeed, interval_centers(int_idx), shipping_factor, depth);
npsd_db_avg = npsd_db_avg + npsd_db_int;
% per区间 maxamp 和 ir_vec
sorted_A_int = preprocessed_data(int_idx).sorted_A;
maxamp_db = max(maxamp_db, 20*log10(max(abs(sorted_A_int))) + sl_db);
ir_vec_int = bharr2ir(preprocessed_data(int_idx).norm_delay, sorted_A_int, Fs); % 列向量
ir_vec_cell{int_idx} = ir_vec_int; % 存 cell
max_ir_len = max(max_ir_len, length(ir_vec_int)); % 更新 max_len
end
npsd_db = npsd_db_avg / num_intervals; %======目前平均也是简化 后续如果可以进行修改=========
nv_db = npsd_db + 10*log10(bw);
% ir_vec 用平均（pad zeros 到 max_len，然后平均）
ir_vec_padded = zeros(max_ir_len, num_intervals) + 1j * zeros(max_ir_len, num_intervals); % 复矩阵
for int_idx = 1:num_intervals
len = length(ir_vec_cell{int_idx});
ir_vec_padded(1:len, int_idx) = ir_vec_cell{int_idx};
end
ir_vec = mean(ir_vec_padded, 2); % 沿列平均，得列向量 (max_len x 1)
snr_dB = uw_isi(ir_vec, sl_db, tx_sig, nv_db, Fs, Fc); % 用平均 ir_vec (列向量)
% 生成着色噪声（不变）
noiseObj = dsp.ColoredNoise('InverseFrequencyPower', 1, ... % α=1 for pink noise，水下低频主导
'SamplesPerFrame', length(tx_sig_channel), ...
'NumChannels', 1, ...
'RandomStream', 'mt19937ar with seed');
colored_noise = noiseObj(); % 生成噪声向量（实数）
colored_noise = colored_noise(:).'; % 强制转置为行向量
% 调整噪声功率匹配nv_db
current_noise_power = var(colored_noise);
target_noise_power = 10^(nv_db / 10);
colored_noise = colored_noise * sqrt(target_noise_power / current_noise_power);
y = tx_sig_channel + colored_noise; % 添加着色噪声
disp(['Strongest tap strength=' num2str(maxamp_db,'%.1f') ' dB; Noise variance=' num2str(nv_db,'%.1f') 'dB']);
disp(['信噪比 = ' num2str(snr_dB,'%.1f') 'dB']);
% === 7. 接收机：OFDM解调 ===
rx_sig = zeros(fftlen, nd); % 原rx_sig（全子载波，调试用）
rx_sig_eq = zeros(num_data_sc, nd); % 均衡后信号（只数据子载波）
% 可选：收集所有符号的H_est和sigma_n2（但per-symbol模式不平均，只用于调试mean）
sigma_n2_all = zeros(1, nd); % 可选收集，用于disp mean
for i = 1:nd
start_idx = (i-1)*sps + 1;
end_idx = i*sps;
symbol_rf = y(start_idx:end_idx); % 提取射频符号
% 下变频到基带
symbol_baseband = symbol_rf .* exp(-1j * 2 * pi * Fc * t_symbol);
% 下采样回基带（从Fs到bw）
symbol_baseband_interp = resample(symbol_baseband, 1, up_factor);
% 确保长度精确
if length(symbol_baseband_interp) > fftlen2
symbol_baseband_interp = symbol_baseband_interp(1:fftlen2);
elseif length(symbol_baseband_interp) < fftlen2
symbol_baseband_interp = [symbol_baseband_interp; zeros(fftlen2 - length(symbol_baseband_interp), 1)];
end
symbol_baseband_interp = symbol_baseband_interp.'; % 转置匹配原代码
% 移除保护间隔
[ich5, qch5] = girem(real(symbol_baseband_interp), imag(symbol_baseband_interp), fftlen2, gilen, 1);
rx = ich5 + 1j * qch5;
% FFT到频域
ry = fft(rx, fftlen) / sqrt(fftlen); % 修改：归一fft
rx_sig(:, i) = ry; % 原存储（可选保留，用于调试）
% --- FDE模块开始（per-symbol） ---
% 步骤1: 提取导频
ry_pilot = ry(pilot_indices);
% % 步骤2: LS信道估计
% H_est_pilot = ry_pilot / pilot_value;
%
% % 步骤3: 插值全H_est
% all_sc = 1:para;
% H_est = interp1(pilot_indices, H_est_pilot, all_sc, 'linear', 'extrap'); % 修改：spline插值
%
% % 加平滑
% H_est = movmean(H_est, 7);
% H_est = H_est*2.5;
% 在接收机FDE循环内，替换步骤2-3（LS+插值）为CS估计
% 步骤2-3: 用OMP稀疏估计（假设时域CIR稀疏）
Phi = exp(-1j * 2 * pi * (0:para-1)' * (0:L-1) / para); % 部分傅里叶字典 (para x L)
Phi_pilot = Phi(pilot_indices, :); % 只取导频行 (num_pilots x L)
% OMP输入：测量值 y_meas = ry_pilot / pilot_value (已知导频除)
y_meas = ry_pilot / pilot_value;
sparsity_level = max([preprocessed_data.num_paths]) * 2; % 稀疏度，2倍最大路径数（裕度）
% 调用OMP（见下方函数）
h_est = omp(y_meas, Phi_pilot, sparsity_level); % h_est: L x 1 稀疏时域CIR
% 转回频域全H_est
H_est = Phi * h_est; % para x 1
% 加平滑（可选，稀疏后轻平滑）
H_est = movmean(H_est, 3); % 窗口3，避免过度平滑稀疏特征
% H_est = H_est;
% 步骤4: 噪声估计（首选null）
ry_null = ry(null_indices);
if use_median_for_noise
sigma_n2 = median(abs(ry_null).^2); % 鲁棒中位数
else
sigma_n2 = mean(abs(ry_null).^2); % 平均
end
% 备选：导频残差（加裕度）
if sigma_n2 < 1e-10
e = ry_pilot - H_est(pilot_indices) * pilot_value;
sigma_n2 = median(abs(e).^2); % 改median
end
sigma_n2 = max(min(sigma_n2 * 1.2, 1e-2), 1e-9); % 加裕度，避免低估
% 步骤5: MMSE均衡（用当前H_est和sigma_n2）
ry_eq = conj(H_est) .* ry ./ (abs(H_est).^2 + sigma_n2);
% ry_eq = ry; %不用均衡
% 加弱子载波丢弃
weak_idx = abs(H_est) < 0.1 * mean(abs(H_est));
ry_eq(weak_idx) = 0; % 修改：丢弃弱载波
% 只取数据子载波
ry_eq_data = ry_eq(data_subcarriers);
rx_sig_eq(:, i) = ry_eq_data;
sigma_n2_all(i) = sigma_n2; % 收集每个符号的噪声方差
% --- FDE模块结束 ---
% 在for i=1:nd循环内，FDE模块后添加
if i == 1
% 计算真实H(f) - 需从preprocessed_data近似（假设用区间中心）
f_sc = (1:para) * subcarrier_spacing + (Fc - bw/2); % 子载波频率轴
H_true = zeros(para,1); % 简化：用最近区间插值
for sc=1:para
[~, int_idx] = min(abs(interval_centers - f_sc(sc)));
norm_delay = preprocessed_data(int_idx).norm_delay;
sorted_A = preprocessed_data(int_idx).sorted_A;
fd_paths = preprocessed_data(int_idx).fd_paths;
for k=1:length(norm_delay)
A_k = sorted_A(k) * 10^(sl_db / 20); % 添加这一行
H_true(sc) = H_true(sc) + A_k * exp(-1j*2*pi*(f_sc(sc)+fd_paths(k))*norm_delay(k));
end
end
figure(15);
plot(abs(H_true), 'b-', 'LineWidth', 1.5); hold on;
plot(abs(H_est), 'r--', 'LineWidth', 1.5);
legend('真实H(f)', '估计H_est');
title('信道幅度响应比较 (第一个符号)');
xlabel('子载波索引'); ylabel('|H|');
end
end
disp(['H_est 平均幅度: ' num2str(mean(abs(H_est)))]);
disp(['H_true 平均幅度: ' num2str(mean(abs(H_true)))]);
% 调试输出平均噪声（可选，使用收集的sigma_n2_all）
disp(['Estimated average sigma_n2 = ' num2str(mean(sigma_n2_all), '%.4f')]);
% QPSK解调（用均衡后信号）
ich6 = real(rx_sig_eq) / kmod;
qch6 = imag(rx_sig_eq) / kmod;
demodata = qpskdemod(ich6, qch6, num_data_sc, nd, ml);
demodata1 = reshape(demodata, 1, num_data_sc*nd*ml);
% === 8. 性能评估 ===
min_len = min(length(source), length(demodata1));
source = source(1:min_len);
signal = demodata1(1:min_len);
[err_number, Pe] = biterr(source, signal);
fprintf('\nOFDM-QPSK性能结果:\n');
fprintf('误码数: %d, 误码率: %f\n', err_number, Pe);
%% 绘图
% 频谱图
figure(6);
N_fft = length(tx_sig);
f = (0:N_fft-1)*(Fs/N_fft);
subplot(3,1,1);
Y_mod = fft(tx_sig, N_fft);
Y_mod_abs = abs(Y_mod/N_fft);
plot(f(1:N_fft/2), 20*log10(Y_mod_abs(1:N_fft/2)));
title('OFDM发射信号幅度谱 (dB)'); xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); grid on;
xlim([0 Fs/2]);
subplot(3,1,2);
Y_tx = fft(tx_sig_channel, N_fft);
Y_tx_abs = abs(Y_tx/N_fft);
plot(f(1:N_fft/2), 20*log10(Y_tx_abs(1:N_fft/2)));
title('OFDM接收信号幅度谱 (dB)'); xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); grid on;
xlim([0 Fs/2]);
sgtitle('OFDM信号频谱对比：发射 vs 接收 (多径+多普勒)');
subplot(3,1,3);
Y_tx = fft(y, N_fft);
Y_tx_abs = abs(Y_tx/N_fft);
plot(f(1:N_fft/2), 20*log10(Y_tx_abs(1:N_fft/2)));
title('OFDM接收信号幅度谱 (dB)'); xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); grid on;
xlim([0 Fs/2]);
sgtitle('OFDM信号频谱对比：发射 vs 接收 (多径+多普勒)');
% 接收端和发射端的信号图绘制
figure(1);
subplot(2,1,1);
stairs(signal);
xlim([1,300]);
ylim([-0.5,1.5]);
title('接收端');
hold on
subplot(2,1,2)
stairs(source);
xlim([1,300]);
ylim([-0.5,1.5]);
title('发送端');
% % 验证着色(粉红)噪声的频谱特性
% figure;
%
% % 生成测试噪声
% test_noise = noiseObj();
%
% % 计算功率谱密度
% [Pxx, F] = pwelch(test_noise, [], [], [], Fs);
%
% % 绘制频谱
% semilogx(F, 10*log10(Pxx));
% xlabel('频率 (Hz)'); ylabel('PSD (dB/Hz)');
% title('粉红噪声频谱 (α=1)');
% grid on;
%
% % 理论参考线（-10 dB/decade = -3 dB/octave）
% hold on;
% theory_slope = -10*log10(F) + 50; % 1/f特性
% semilogx(F, theory_slope, 'r--', 'LineWidth', 1.5);
% legend('实际噪声', '理论1/f斜率', 'Location', 'best');
%%
%% 函数
% 2PSK调制
function [ich, qch] = qpskmod(paradata, para, nd, ml)
ich = zeros(para, nd);
qch = zeros(para, nd);
for i = 1:para
for j = 1:nd
bits = paradata(i, (j-1)*ml+1:j*ml);
if isequal(bits, [0 0])
ich(i,j) = 1; qch(i,j) = 1;
elseif isequal(bits, [0 1])
ich(i,j) = 1; qch(i,j) = -1;
elseif isequal(bits, [1 1])
ich(i,j) = -1; qch(i,j) = -1;
elseif isequal(bits, [1 0])
ich(i,j) = -1; qch(i,j) = 1;
end
end
end
end
%2PSK解调
function demodata = qpskdemod(ich, qch, para, nd, ml)
demodata = zeros(para, nd*ml);
for i = 1:para
for j = 1:nd
if ich(i,j) > 0 && qch(i,j) > 0
demodata(i, (j-1)*ml+1:j*ml) = [0 0];
elseif ich(i,j) > 0 && qch(i,j) < 0
demodata(i, (j-1)*ml+1:j*ml) = [0 1];
elseif ich(i,j) < 0 && qch(i,j) < 0
demodata(i, (j-1)*ml+1:j*ml) = [1 1];
elseif ich(i,j) < 0 && qch(i,j) > 0
demodata(i, (j-1)*ml+1:j*ml) = [1 0];
end
end
end
end
%添加保护间隔 循环前缀
function [ich3, qch3] = giins(ich2, qch2, fftlen, gilen, nd)
ich3 = zeros(fftlen + gilen, nd);
qch3 = zeros(fftlen + gilen, nd);
for i = 1:nd
ich3(1:gilen, i) = ich2(end-gilen+1:end, i); % 循环前缀
ich3(gilen+1:end, i) = ich2(:, i);
qch3(1:gilen, i) = qch2(end-gilen+1:end, i);
qch3(gilen+1:end, i) = qch2(:, i);
end
end
%保护间隔移除
function [ich5, qch5] = girem(ich4, qch4, fftlen2, gilen, nd)
ich5 = ich4(gilen+1:end, :);
qch5 = qch4(gilen+1:end, :);
end
% 转未归一化离散时间冲击响应
function ir_vec = bharr2ir(norm_delay, sorted_A, Fs)
delay_max = max(norm_delay);
cir_length = round(delay_max*Fs);
ir_vec = zeros(cir_length+1, 1) + 1j * zeros(cir_length + 1, 1);
for icn = 1: length(norm_delay)
arr_id_no = round(norm_delay(icn)*Fs);
if arr_id_no >= 0
ir_vec(arr_id_no + 1) = ir_vec(arr_id_no + 1) + sorted_A(icn);
end
end
end
% 海洋背景噪声函数
function [npsd_db] = improved_ambientnoise_psd_v2(windspeed, Fc, shipping_factor, depth)
if nargin < 3, shipping_factor = 0.5; end
if nargin < 4, depth = 0; end
f = Fc / 1000; % f in kHz
% 湍流噪声 (体积源，无深度衰减)
ANturb_dB = 17 - 30 * log10(f);
% 风浪噪声 (表面源，深度衰减)
ws = min(windspeed, 20);
ANwind_dB = 50 + 7.5 * sqrt(ws) + 20 * log10(f) - 40 * log10(f + 0.4);
if depth > 0
surface_att = 0.01 * depth * f; % 经验衰减 (dB/m * kHz scaling)
ANwind_dB = ANwind_dB - surface_att;
end
% 热噪声 (体积源，无深度衰减)
ANthermo_dB = -15 + 20 * log10(f);
% 船舶噪声 (表面源，深度衰减)
ANship_dB = 40 + 20 * (shipping_factor - 0.5) + 26 * log10(f) - 60 * log10(f + 0.03);
if depth > 0
surface_att = 0.01 * depth * f; % 同上
ANship_dB = ANship_dB - surface_att;
end
% 总线性功率
total_power = 10^(ANturb_dB/10) + 10^(ANwind_dB/10) + 10^(ANthermo_dB/10) + 10^(ANship_dB/10);
npsd_db = 10 * log10(total_power); % dB re 1 μPa²/Hz
end
% 求真实信噪比
function snr_dB = uw_isi(ir_vec_no, sl_db, sig_modulate, nv_db, Fs, Fc)
t_ir_no = ((0:length(ir_vec_no) - 1) / Fs)'; % 时间向量（列）
adj_ir_complex = ir_vec_no * 10^(sl_db / 20); % 绝对幅度，复基带
rf_adj_complex = adj_ir_complex .* exp(1j * 2 * pi * Fc * t_ir_no); % 上变频
adj_ir = real(rf_adj_complex); % 射频实 CIR
tx_sig = conv(sig_modulate, adj_ir); % 实 conv
tx_sig = tx_sig(1:length(sig_modulate)); % 对齐
tx_sig = reshape(tx_sig, 1, []); % 新增：强制为 1×L 行向量
namp = 10^(nv_db/20); % 噪声的强度值
noise = namp * randn(size(tx_sig)); % 实高斯噪声
snr_dB = 10*log10(mean(tx_sig.^2)/mean(noise.^2)); % 平均功率信噪比
end
% 在代码末尾添加OMP函数（简单实现，无需工具箱）
function x = omp(y, A, K)
% OMP: Orthogonal Matching Pursuit
% y: measurements (m x 1)
% A: dictionary (m x n)
% K: sparsity level
[m, n] = size(A);
x = zeros(n, 1);
r = y; % 残差
supp = []; % 支持集
for k = 1:K
corr = abs(A' * r); % 相关度
[~, idx] = max(corr); % 最相关原子
supp = [supp, idx]; % 添加支持
As = A(:, supp); % 子字典
xs = As \ y; % LS求解
r = y - As * xs; % 更新残差
if norm(r) < 1e-6, break; end % 早停
end
x(supp) = xs;
end
理解这部分代码 其中preprocessed_data为：[5.278322e-06 + 0.000000i,4.768421e-05 + 0.000000i,1.475299e-05 + 2.252722e-05i,-4.525621e-05 - 3.956423e-12i] [0,1.9073486e-06,0.23500347,0.41591835] [8.3433142,8.3356848,17.706266,-15.234738] [-3.1621485,-3.1685233,-15.966362,13.235704] [0;0;0;0] 4
[6.152235e-06 + 0.000000i,4.680981e-05 + 0.000000i,2.764621e-06 + 4.217828e-06i,1.184823e-05 + 1.800473e-05i,1.835814e-07 + 2.778731e-07i,-4.525598e-05 - 3.956404e-12i] [0,1.9073486e-06,0.23499966,0.23500156,0.23500347,0.41591930] [8.3428640,8.3355999,17.702553,17.707031,17.712931,-15.234738] [-3.1625235,-3.1685936,-15.957605,-15.968169,-15.982075,13.235705] [0;0;0;0;0;0] 6
[3.738785e-07 + 0.000000i,5.258722e-05 + 0.000000i,1.475749e-05 + 2.252527e-05i,-1.245270e-05 - 1.088649e-12i,-3.280272e-05 - 2.867705e-12i] [0,2.8610229e-06,0.23500347,0.41591930,0.41592216] [8.3467369,8.3363705,17.706266,-15.231303,-15.236040] [-3.1592860,-3.1679497,-15.966364,13.235600,13.235744] [0;0;0;0;0] 5
[5.296488e-05 + 0.000000i,8.958914e-06 + 1.366344e-05i,5.827515e-06 + 8.843062e-06i,-2.713330e-07 - 2.372068e-14i,-4.498314e-05 - 3.932551e-12i] [0,0.23500156,0.23500443,0.41591835,0.41592121] [8.3364449,17.704702,17.708681,-15.226946,-15.234783] [-3.1678879,-15.962675,-15.972057,13.235473,13.235704] [0;0;0;0;0] 5
[9.376385e-06 + 0.000000i,4.358450e-05 + 0.000000i,1.657540e-07 + 2.529105e-07i,1.462083e-05 + 2.225388e-05i,-2.073672e-06 - 1.812862e-13i,-4.317944e-05 - 3.774866e-12i] [0,1.9073486e-06,0.23499870,0.23500252,0.41591644,0.41592026] [8.3416624,8.3353233,17.699976,17.706341,-15.229079,-15.235007] [-3.1635282,-3.1688249,-15.951522,-15.966533,13.235535,13.235712] [0;0;0;0;0;0] 6
[1.130504e-07 + 0.000000i,5.284886e-05 + 0.000000i,1.476881e-05 + 2.251722e-05i,-4.525124e-05 - 3.955989e-12i] [0,2.8610229e-06,0.23500633,0.41591835] [8.3473539,8.3364229,17.706270,-15.234735] [-3.1587691,-3.1679072,-15.966367,13.235701] [0;0;0;0] 4
[1.202769e-05 + 0.000000i,4.093478e-05 + 0.000000i,1.914058e-06 + 2.919557e-06i,1.286570e-05 + 1.955519e-05i,1.697104e-08 + 2.567433e-08i,-4.525042e-05 - 3.955918e-12i] [0,1.9073486e-06,0.23499966,0.23500252,0.23500443,0.41591930] [8.3409615,8.3351183,17.702517,17.706816,17.713816,-15.234735] [-3.1641140,-3.1689966,-15.957516,-15.967661,-15.984159,13.235702] [0;0;0;0;0;0] 6
[5.296262e-05 + 0.000000i,1.476896e-05 + 2.251103e-05i,2.184308e-09 + 3.302675e-09i,-4.525127e-05 - 3.955991e-12i] [0,0.23500252,0.23500443,0.41591740] [8.3364439,17.706268,17.714851,-15.234735] [-3.1678874,-15.966363,-15.986600,13.235703] [0;0;0;0] 4
[2.531477e-07 + 0.000000i,5.270921e-05 + 0.000000i,1.268601e-05 + 1.934451e-05i,2.090223e-06 + 3.169412e-06i,-4.525162e-05 - 3.956023e-12i] [0,1.9073486e-06,0.23500156,0.23500443,0.41591930] [8.3458652,8.3363981,17.705692,17.709766,-15.234737] [-3.1600146,-3.1679258,-15.965010,-15.974618,13.235703] [0;0;0;0;0] 5
[1.669531e-05 + 0.000000i,3.626637e-05 + 0.000000i,1.476279e-05 + 2.252063e-05i,-4.525160e-05 - 3.956020e-12i] [0,2.8610229e-06,0.23500347,0.41591740] [8.3400440,8.3347864,17.706266,-15.234734] [-3.1648798,-3.1692729,-15.966362,13.235701] [0;0;0;0] 4
[5.296050e-05 + 0.000000i,1.303849e-05 + 1.987655e-05i,1.738836e-06 + 2.636477e-06i,-4.525128e-05 - 3.955992e-12i] [0,0.23500061,0.23500156,0.41591740] [8.3364439,17.705793,17.709822,-15.234735] [-3.1678884,-15.965250,-15.974747,13.235702] [0;0;0;0] 4
[5.242816e-05 + 0.000000i,5.367077e-07 + 0.000000i,5.511806e-06 + 8.405864e-06i,9.245658e-06 + 1.404052e-05i,3.630870e-08 + 5.496389e-08i,-1.489615e-05 - 1.302262e-12i,-4.853558e-07 - 4.243115e-14i,-2.986919e-05 - 2.611248e-12i] [0,0,0.23499966,0.23500061,0.23500347,0.41591549,0.41591740,0.41591930] [8.3363619,8.3445988,17.704155,17.707502,17.712618,-15.237297,-15.228613,-15.233560] [-3.1679575,-3.1610742,-15.961385,-15.969279,-15.981338,13.235780,13.235521,13.235667] [0;0;0;0;0;0;0;0] 8
[6.777499e-07 + 0.000000i,2.132690e-05 + 0.000000i,3.095676e-05 + 0.000000i,8.880509e-06 + 1.353449e-05i,5.904511e-06 + 8.961507e-06i,5.309370e-09 + 8.033200e-09i,-2.889215e-05 - 2.525832e-12i,-1.635774e-05 - 1.430039e-12i] [0,1.9073486e-06,2.8610229e-06,0.23500156,0.23500347,0.23500538,0.41591740,0.41592121] [8.3442001,8.3391142,8.3344374,17.704956,17.708237,17.713577,-15.233394,-15.237106] [-3.1614068,-3.1656585,-3.1695650,-15.963274,-15.971011,-15.983598,13.235660,13.235777] [0;0;0;0;0;0;0;0] 8
[8.497629e-07 + 0.000000i,5.211246e-05 + 0.000000i,6.154986e-06 + 9.384284e-06i,8.638848e-06 + 1.311767e-05i,-4.525769e-05 - 3.956553e-12i] [0,1.9073486e-06,0.23500156,0.23500252,0.41591930] [8.3438177,8.3363247,17.704391,17.707605,-15.234735] [-3.1617262,-3.1679881,-15.961941,-15.969522,13.235703] [0;0;0;0;0] 5
[5.296258e-05 + 0.000000i,6.295411e-07 + 9.592380e-07i,1.416199e-05 + 2.153556e-05i,3.521108e-09 + 5.327710e-09i,-4.525635e-05 - 3.956436e-12i] [0,0.23500061,0.23500347,0.23500538,0.41591930] [8.3364439,17.702106,17.706450,17.713511,-15.234736] [-3.1678872,-15.956551,-15.966798,-15.983438,13.235703] [0;0;0;0;0] 5
[2.761986e-05 + 0.000000i,2.534269e-05 + 0.000000i,1.475420e-05 + 2.248970e-05i,1.779290e-08 + 2.693569e-08i,-4.525455e-05 - 3.956278e-12i] [0,4.7683716e-06,0.23500538,0.23500729,0.41591930] [8.3386183,8.3340750,17.706263,17.712555,-15.234735] [-3.1660714,-3.1698675,-15.966346,-15.981189,13.235702] [0;0;0;0;0] 5
[2.842072e-07 + 0.000000i,5.255642e-05 + 0.000000i,1.214780e-07 + 0.000000i,1.477458e-05 + 2.250728e-05i,2.308047e-09 + 3.492371e-09i,-4.525803e-05 - 3.956582e-12i] [0,3.8146973e-06,4.7683716e-06,0.23500443,0.23500633,0.41592026] [8.3444157,8.3364239,8.3277340,17.706266,17.713448,-15.234736] [-3.1612265,-3.1679077,-3.1751609,-15.966364,-15.983290,13.235702] [0;0;0;0;0;0] 6
[5.100477e-05 + 0.000000i,1.956443e-06 + 0.000000i,7.479953e-06 + 1.139905e-05i,7.301478e-06 + 1.108490e-05i,1.223933e-08 + 1.852872e-08i,-4.525305e-05 - 3.956147e-12i] [0,2.8610229e-06,0.23500252,0.23500443,0.23500633,0.41592121] [8.3362131,8.3424339,17.704800,17.707764,17.712526,-15.234738] [-3.1680801,-3.1628835,-15.962905,-15.969894,-15.981122,13.235703] [0;0;0;0;0;0] 6
[4.708412e-07 + 0.000000i,5.249402e-05 + 0.000000i,1.476952e-05 + 2.251604e-05i,-4.422494e-05 - 3.866267e-12i,-1.028481e-06 - 8.991268e-14i] [0,9.5367432e-07,0.23500252,0.41591644,0.41591835] [8.3437347,8.3363800,17.706264,-15.234846,-15.229898] [-3.1617975,-3.1679425,-15.966363,13.235704,13.235559] [0;0;0;0;0] 5
[3.491842e-05 + 0.000000i,1.804234e-05 + 0.000000i,1.476812e-05 + 2.250106e-05i,8.321412e-09 + 1.259769e-08i,-4.525337e-05 - 3.956176e-12i] [0,3.8146973e-06,0.23500347,0.23500633,0.41591930] [8.3379202,8.3335915,17.706264,17.712500,-15.234734] [-3.1666546,-3.1702724,-15.966356,-15.981060,13.235702] [0;0;0;0;0] 5
你检查一下这个通讯模拟有什么问题 注释中标注不准的忽视 多普勒频移为0的问题忽略 重点关注我的OFMD的参数设置和用以补偿多径效应的信道估计和均衡器 同时评估一下 我现在的这个均衡器能否有效地解决当前多径效应严重的问题
最新发布
扫一扫
2900
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
