% 读取实部和虚部COE文件并构建复数信号的完整代码
% 文件目录: G:\matlab.2\work
% 说明：echo_dbf1_s_i = 实部数据，echo_dbf1_s_q = 虚部数据，最终复数信号 = echo_dbf1_s
% ------------------------------
% 第一步：读取实部COE文件（echo_dbf1_s_i.coe）
% ------------------------------
file_path = 'G:\matlab.2\work\';
filename_real = [file_path, 'echo_dbf1_s_i.coe'];
echo_dbf1_s_i = readCoeFile(filename_real);
disp('实部COE文件（echo_dbf1_s_i.coe）读取完成！');
fprintf('实部数据总数: %d\n', length(echo_dbf1_s_i));
disp('实部前20个数据:');
for i = 1:min(20, length(echo_dbf1_s_i))
fprintf('echo_dbf1_s_i[%d] = %d\n', i-1, echo_dbf1_s_i(i));
end
figure('Name', '实部数据时域图（echo_dbf1_s_i）');
subplot(2,1,1);
plot(echo_dbf1_s_i);
title('实部COE文件数据（echo_dbf1_s_i）- 时域图');
xlabel('样本索引');
ylabel('实部幅值');
grid on;
subplot(2,1,2);
stem(echo_dbf1_s_i(1:min(100, length(echo_dbf1_s_i))), 'filled');
title('实部COE文件数据（echo_dbf1_s_i）- 前100个样本');
xlabel('样本索引');
ylabel('实部幅值');
grid on;
% ------------------------------
% 第二步：读取虚部COE文件（echo_dbf1_s_q.coe）
% ------------------------------
filename_imag = [file_path, 'echo_dbf1_s_q.coe'];
echo_dbf1_s_q = readCoeFile(filename_imag);
disp('虚部COE文件（echo_dbf1_s_q.coe）读取完成！');
fprintf('虚部数据总数: %d\n', length(echo_dbf1_s_q));
disp('虚部前20个数据:');
for i = 1:min(20, length(echo_dbf1_s_q))
fprintf('echo_dbf1_s_q[%d] = %d\n', i-1, echo_dbf1_s_q(i));
end
figure('Name', '虚部数据时域图（echo_dbf1_s_q）');
subplot(2,1,1);
plot(echo_dbf1_s_q);
title('虚部COE文件数据（echo_dbf1_s_q）- 时域图');
xlabel('样本索引');
ylabel('虚部幅值');
grid on;
subplot(2,1,2);
stem(echo_dbf1_s_q(1:min(100, length(echo_dbf1_s_q))), 'filled');
title('虚部COE文件数据（echo_dbf1_s_q）- 前100个样本');
xlabel('样本索引');
ylabel('虚部幅值');
grid on;
% ------------------------------
% 第三步：构建复数信号
% ------------------------------
if length(echo_dbf1_s_i) ~= length(echo_dbf1_s_q)
error('实部（echo_dbf1_s_i）和虚部（echo_dbf1_s_q）数据长度不一致，请检查COE文件');
end
echo_dbf1_s = echo_dbf1_s_i + 1j * echo_dbf1_s_q;
% ------------------------------
% 第四步：复数信号波形展示
% ------------------------------
figure('Name', '复数信号（echo_dbf1_s）时域波形', 'Position', [100, 200, 1000, 600]);
% 4.1 完整时域波形
subplot(2,1,1);
plot(1:length(echo_dbf1_s), real(echo_dbf1_s), 'b-', 'LineWidth', 1.2, 'DisplayName', '实部');
hold on;
plot(1:length(echo_dbf1_s), imag(echo_dbf1_s), 'r--', 'LineWidth', 1.2, 'DisplayName', '虚部');
hold off;
title('复数信号（echo_dbf1_s）完整时域波形');
xlabel('样本索引');
ylabel('幅值');
grid on;
legend();
% ------------------------------
% 第五步：STFT 分析、镂空处理及逆 STFT 变换
% ------------------------------
% 参数设置
c = 3e8; % 光速（备用参数）
fs_base = 250e6; % 采样频率
win = hamming(128); % 汉宁窗（长度32）
overlap = 64 % 帧重叠长度
Nfft = 128; % FFT点数
STFT_gate = 220; % 镂空处理阈值
% 定义信号段数（若为单段信号，N_trans=1）
N_trans = 1;
% 初始化存储变量（维度与信号匹配）
signal_len = length(echo_dbf1_s);
sss_sum_u = zeros(N_trans, signal_len); % 原始信号存储
sss_stft_sum_u = zeros(N_trans, signal_len); % 逆STFT结果（s_sum_u）
% ------------------------------
% 3. 执行STFT变换（手动实现，替换内置stft）
% ------------------------------
wlen = length(win); % 128
hop = overlap; % 64
nfft = Nfft; % 128
h = win; % 汉明窗
x1 = echo_dbf1_s(:); % 转为列向量
signal_length = length(x1);
num_frames = floor((signal_length - wlen) / hop) + 1;
% 初始化STFT结果矩阵 (nfft × num_frames)
s_sum_u = zeros(nfft, num_frames);
% 手动实现STFT（支持复数）
for frame_idx = 1:num_frames
start_idx = (frame_idx - 1) * hop + 1;
end_idx = start_idx + wlen - 1;
if end_idx <= signal_length
frame_data = x1(start_idx:end_idx);
windowed_data = frame_data .* h; % 复数 × 实窗
s_sum_u(:, frame_idx) = fft(windowed_data, nfft);
end
end
% 注意：MATLAB stft 返回的是 'frequency × time'，且默认 fftshift 未应用
% 所以我们不需要 fftshift
% 计算时间轴 t（秒）
t = (0:num_frames-1) * hop / fs_base;
% 计算频率轴 f（Hz）—— 与 MATLAB stft 对齐（未 fftshift，从 0 到 fs）
% MATLAB 的 stft 默认频率范围是 [0, fs)，共 nfft 点
f = (0:nfft-1) * fs_base / nfft;
% 5. 镂空处理：超过阈值的频点置0
for i = 1:63
if abs(s_sum_u(1, i)) > STFT_gate
s_sum_u(1,i) = 0;
end
end
% 6. 逆STFT变换（手动实现，替换内置istft）
% ========================
% 手动 ISTFT 代码开始
% ========================
% 6.1 准备输出信号
output_length = (num_frames - 1) * hop + wlen;
x_recon = zeros(output_length, 1); % 重建信号（叠加结果）
window_sum = zeros(output_length, 1); % 窗函数叠加权重（用于归一化）
% 6.2 逐帧重建
for frame_idx = 1:num_frames
% 1. 从频域做 IFFT 得到时域帧（复数）
frame_spectrum = s_sum_u(:, frame_idx);
frame_time = ifft(frame_spectrum, nfft); % 注意：这里用了 ifft，但没用 istft
frame_time = frame_time(1:wlen); % 只取窗长部分（假设 nfft >= wlen）
% 2. 计算该帧在输出信号中的起始/结束位置
start_idx = (frame_idx - 1) * hop + 1;
end_idx = start_idx + wlen - 1;
% 3. 重叠相加（OLA: OverLap-Add）
if end_idx <= output_length
x_recon(start_idx:end_idx) = x_recon(start_idx:end_idx) + frame_time;
window_sum(start_idx:end_idx) = window_sum(start_idx:end_idx) + h;
end
end
% 6.3 归一化（避免窗函数重叠导致的幅度失真）
% 注意：理想情况下，窗函数在 hop 步长下应满足 "常数重叠相加"（COLA）性质
valid_idx = window_sum > 1e-10; % 避免除零
x_recon(valid_idx) = x_recon(valid_idx) ./ window_sum(valid_idx);
% 最终重建信号（实部，因为原始信号是实或复，但STFT是复数）
istft_result = x_recon; % 若原始为复信号，保留复数；若为实信号，可取 real(x_recon)
% ========================
% 手动 ISTFT 代码结束
% ========================
% ------------------------------
% 重新绘制镂空后的 STFT 信号波形图（第1频率 bin）
% 横坐标：对应的原始信号采样点序号
% 纵坐标：STFT 系数的实部 / 虚部
% ------------------------------
% 提取镂空处理后的 STFT 第1个频率 bin（已修改过的 s_sum_u）
stft_bin_processed = s_sum_u(1, :); % 复数向量，长度 = 时间帧数
% 时间轴 t 由 stft 函数返回（单位：秒）
% 转换为对应的原始信号采样点位置（取每帧中心点）
sample_points = t * fs_base + 1; % 转为1起始索引
sample_points = round(sample_points); % 取整为采样点序号
% 确保长度一致（防止边界误差）
N_plot = min(length(stft_bin_processed), length(sample_points));
stft_bin_processed = stft_bin_processed(1:N_plot);
sample_points = sample_points(1:N_plot);
% 绘制实部和虚部
figure('Name', '镂空后 STFT 系数波形图（第1频率 bin）', 'Position', [300, 120, 1000, 600]);
subplot(2,1,1);
plot(sample_points, real(stft_bin_processed), 'b-o', 'MarkerFaceColor', 'b', 'MarkerSize', 4);
title('镂空后 STFT 系数的实部（第1频率 bin）');
xlabel('采样点序号');
ylabel('实部');
grid on;
subplot(2,1,2);
plot(sample_points, imag(stft_bin_processed), 'r-o', 'MarkerFaceColor', 'r', 'MarkerSize', 4);
title('镂空后 STFT 系数的虚部（第1频率 bin）');
xlabel('采样点序号');
ylabel('虚部');
grid on;
% ------------------------------
% 第六步：绘制 istft 重建信号的实部和虚部
% ------------------------------
figure('Name', 'ISTFT 重建信号的实部与虚部', 'Position', [200, 150, 1000, 600]);
% 绘制实部
subplot(2, 1, 1);
plot(1:length(istft_result), real(istft_result), 'b-', 'LineWidth', 1.2);
title('ISTFT 重建信号 - 实部');
xlabel('样本索引');
ylabel('实部幅值');
grid on;
% 绘制虚部
subplot(2, 1, 2);
plot(1:length(istft_result), imag(istft_result), 'r--', 'LineWidth', 1.2);
title('ISTFT 重建信号 - 虚部');
xlabel('样本索引');
ylabel('虚部幅值');
grid on;
% ------------------------------
% COE文件读取函数
% ------------------------------
function data = readCoeFile(filename)
% 读取Xilinx标准COE文件（支持2/10/16进制）
% 输入：filename - COE文件完整路径
% 输出：data - 十进制数据向量
if ~exist(filename, 'file')
error('文件不存在：%s', filename);
end
fprintf('正在读取文件：%s\n', filename);
fid = fopen(filename, 'r');
if fid == -1
error('无法打开文件（权限或路径错误）：%s', filename);
end
try
radix = 10; % 默认十进制
data = [];
reading_data = false;
line_num = 0;
while ~feof(fid)
line = fgetl(fid);
line_num = line_num + 1;
if isempty(line), continue; end
line = strtrim(line);
if startsWith(line, ';'), continue; end % 跳过注释
line_lower = lower(line);
% 解析数据基数
if ~reading_data && contains(line_lower, 'memory_initialization_radix')
radix_match = regexp(line_lower, '=\s*(\d+)', 'tokens');
if ~isempty(radix_match)
radix = str2double(radix_match{1}{1});
if ~ismember(radix, [2, 10, 16])
error('第%d行：不支持的基数%d（仅支持2/10/16进制）', line_num, radix);
end
fprintf('检测到数据基数：%d进制\n', radix);
end
continue;
end
% 解析数据起始标记
if contains(line_lower, 'memory_initialization_vector')
reading_data = true;
fprintf('进入数据读取阶段...\n');
if contains(line, '=')
data_part = extractAfter(line, '=');
data_part = strtrim(data_part);
if endsWith(data_part, ';'), data_part = data_part(1:end-1); end
new_data = parseDataLine(data_part, radix);
data = [data; new_data];
end
continue;
end
% 读取数据行
if reading_data
if endsWith(line, ';'), line = line(1:end-1); end
new_data = parseDataLine(line, radix);
data = [data; new_data];
end
end
if isempty(data)
error('未从文件中读取到有效数据，请检查COE格式');
end
catch err
error('读取失败（第%d行）：%s', line_num, err.message);
finally
fclose(fid);
end
end
% ------------------------------
% 辅助函数：解析单行COE数据
% ------------------------------
function data_line = parseDataLine(line, radix)
data_str = strsplit(line, ',');
data_str = strtrim(data_str);
data_str = data_str(~cellfun('isempty', data_str));
data_line = zeros(length(data_str), 1);
for i = 1:length(data_str)
item = data_str{i};
switch radix
case 2, data_line(i) = bin2dec(item);
case 10, data_line(i) = str2double(item);
case 16, data_line(i) = hex2dec(upper(item));
end
end
end分析代码为什么，镂空后，istft输出一个0都没有
最新发布
本文介绍了一个简单的C语言程序，用于计算从1到指定整数m的所有整数的平方和。通过for循环逐个计算每个整数的平方并累加到变量sun中，最后输出结果。
扫一扫
1130
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
