论坛 › 机器人GPT › ZBF

Anonymous 发表于 2025-1-6 09:36:16

还是要先对，zbf读出的初始的数据进行归一化处理，要不然读出的fft数据还是会有一堆NaN

机器人Gpt 发表于 2025-1-6 09:36:30

LIYU 发表于 2025-01-06 09:36
还是要先对，zbf读出的初始的数据进行归一化处理，要 ...

为了避免在傅里叶变换过程中出现 `NaN` 或 `Inf`，对初始数据进行适当的预处理是很有必要的。我们可以在光场构建之前进行归一化处理，并清除可能的异常值。以下是如何执行这些步骤的完整代码示例，确保所有数据的数值范围合理：

% 打开 ZBF 文件
filename = 'your_file.zbf'; % 请替换为实际的文件名
fileID = fopen(filename, 'r', 'l'); % 'l' denotes little-endian

if fileID == -1
    error('无法打开文件');
end

% 假设我们知道数据大小为 Nx x Ny
Nx = 256; % 从文件中解析
Ny = 256; % 从文件中解析
dataSize = Nx * Ny;

% 读取照度数据
intensity = fread(fileID, dataSize, 'double');

% 读取相位数据
phase = fread(fileID, dataSize, 'double');

fclose(fileID);

% 对照度和相位数据进行预处理
% 将 NaN 和 Inf 替换为有限值
intensity(isnan(intensity) | isinf(intensity)) = 0;
phase(isnan(phase) | isinf(phase)) = 0;

% 对照度数据进行归一化处理
intensity = intensity / max(intensity(:)); % 归一化至

% 确保数据组织成矩阵形式
intensity = reshape(intensity, Nx, Ny);
phase = reshape(phase, Nx, Ny);

% 形成复数形式的光场
complexField = sqrt(intensity) .* exp(1i * phase);

% 修正 potential NaNs/Infs in complexField after combination
complexField(isnan(complexField) | isinf(complexField)) = 0;

% 傅里叶变换
fftResult = fftshift(fft2(complexField));

% 能量分布
energyDistribution = abs(fftResult).^2;

% 计算中心坐标
centerX = ceil(Nx / 2);
centerY = ceil(Ny / 2);

% 网格坐标
= meshgrid(1:Nx, 1:Ny);

% 计算径向平均
maxRadius = min(centerX, centerY);
radialProfile = zeros(maxRadius, 1);
count = zeros(maxRadius, 1);

% 收集每个半径上的能量值
for x = 1:Nx
    for y = 1:Ny
      radius = round(sqrt((x - centerX)^2 + (y - centerY)^2));
      if radius >= 1 && radius <= maxRadius
            radialProfile(radius) = radialProfile(radius) + energyDistribution(x, y);
            count(radius) = count(radius) + 1;
      end
    end
end

% 归一化径向能量分布
radialProfile = radialProfile ./ count;
radialProfile = radialProfile / max(radialProfile);

% 绘制径向能量分布图
figure;
plot(1:maxRadius, radialProfile);
xlabel('Radius in Frequency Domain');
ylabel('Normalized Energy');
title('Radial Energy Distribution');

### 说明：

1. **预处理步骤**：在读取后立即对原始数据进行归一化，并处理任何无效值（如 `NaN` 或 `Inf`），这有助于确保傅里叶变换输入数据的稳定性。

2. **径向平均**：计算径向能量分布并绘制，以揭示频域中的对称性和能量特征。

3. **验证和调试**：确保所有中间数据和步骤都经过验证，没有异常值，这将有助于分析和调试过程。

通过这种方式，你应该可以获得一个更为稳定和可解释的频域能量分布结果。
--- 光学专家Gpt

Anonymous 发表于 2025-1-6 10:18:50

绘图还是不对啊，我需要得到与圆形能量分布图的差异：而圆形的能量分布图通常是在极坐标下绘制，它能够直观地展示在不同角度和半径下的能量分布情况。当前代码没有体现角度相关的信息，所以不能得到圆形的能量分布图。

机器人Gpt 发表于 2025-1-6 10:19:00

LIYU 发表于 2025-01-06 10:18
绘图还是不对啊，我需要得到与圆形能量分布图的差异： ...

理解你的需求后，若要获得一个在二维频域中能量分布的图，并体现出各方向的对称性（或缺乏对称性），我们可以通过绘制二维能量谱图来实现。这种图像通常通过展示能够直观反映出傅里叶变换后信号的频率结构。

这里是如何实现这一点的代码：

% 打开 ZBF 文件
filename = 'your_file.zbf'; % 请替换为实际的文件名
fileID = fopen(filename, 'r', 'l'); % 'l' denotes little-endian

if fileID == -1
    error('无法打开文件');
end

% 假设我们知道数据大小为 Nx x Ny
Nx = 256; % 从文件中解析
Ny = 256; % 从文件中解析
dataSize = Nx * Ny;

% 读取照度数据
intensity = fread(fileID, dataSize, 'double');

% 读取相位数据
phase = fread(fileID, dataSize, 'double');

fclose(fileID);

% 对照度和相位数据进行预处理
intensity(isnan(intensity) | isinf(intensity)) = 0;
phase(isnan(phase) | isinf(phase)) = 0;

% 对照度数据进行归一化处理
intensity = intensity / max(intensity(:));

% 确保数据组织成矩阵形式
intensity = reshape(intensity, Nx, Ny);
phase = reshape(phase, Nx, Ny);

% 形成复数形式的光场
complexField = sqrt(intensity) .* exp(1i * phase);

% 修正可能的 NaNs/Infs 后组合
complexField(isnan(complexField) | isinf(complexField)) = 0;

% 傅里叶变换
fftResult = fftshift(fft2(complexField));

% 能量分布
energyDistribution = abs(fftResult).^2;

% 归一化能量分布
normalizedEnergyDistribution = energyDistribution / max(energyDistribution(:));

% 使用 imagesc 绘制二维能量分布图
figure;
imagesc(normalizedEnergyDistribution);
colormap('jet');
colorbar;
title('2D Frequency Domain Energy Distribution');
xlabel('Spatial Frequency X');
ylabel('Spatial Frequency Y');
axis image;

### 说明：
- **二维频域图**：`imagesc` 函数用于直接可视化二维能量分布，横轴和纵轴代表空间频率的两个方向。

- **对称性观察**：通过这个图，你可以观察能量在频域中的分布样式，以及是否存在圆形对称或其他形式的对称性。

- **数据处理**：在傅里叶变换之前已经进行了数据预处理，以消除任何 `NaN` 和 `Inf` 的影响，同时进行了强度归一化以避免数值问题。

这种方法让你不需要转换到极坐标系，就能从直角坐标系下查看能量的分布特征。能量分布的形状将指示输入信号的频率特性，特别是其对称性或非对称性。
--- 光学专家Gpt

Anonymous 发表于 2025-1-6 10:22:48

但是这种方法我没办法直观地从角度来看频数，再得到对应角度的能量，占据总能量的百分比，能修改一下吗

机器人Gpt 发表于 2025-1-6 10:23:00

LIYU 发表于 2025-01-06 10:22
但是这种方法我没办法直观地从角度来看频数，再得到对 ...

要得到频域能量在不同角度下的分布，并计算其占据总能量的百分比，可以将傅里叶变换的结果转换为极坐标系，然后进行径向和角度上的积分。以下是如何实现这些步骤的代码：

% 打开 ZBF 文件
filename = 'your_file.zbf'; % 请替换为实际的文件名
fileID = fopen(filename, 'r', 'l'); % 'l' denotes little-endian

if fileID == -1
    error('无法打开文件');
end

% 假设我们知道数据大小为 Nx x Ny
Nx = 256; % 从文件中解析
Ny = 256; % 从文件中解析
dataSize = Nx * Ny;

% 读取照度数据
intensity = fread(fileID, dataSize, 'double');

% 读取相位数据
phase = fread(fileID, dataSize, 'double');

fclose(fileID);

% 对照度和相位数据进行预处理
intensity(isnan(intensity) | isinf(intensity)) = 0;
phase(isnan(phase) | isinf(phase)) = 0;

% 对照度数据进行归一化处理
intensity = intensity / max(intensity(:));

% 确保数据组织成矩阵形式
intensity = reshape(intensity, Nx, Ny);
phase = reshape(phase, Nx, Ny);

% 形成复数形式的光场
complexField = sqrt(intensity) .* exp(1i * phase);

% 修正可能的 NaNs/Infs 后组合
complexField(isnan(complexField) | isinf(complexField)) = 0;

% 傅里叶变换
fftResult = fftshift(fft2(complexField));

% 能量分布
energyDistribution = abs(fftResult).^2;
totalEnergy = sum(energyDistribution(:));

% 创建网格来计算半径和角度
= meshgrid(1:Nx, 1:Ny);
centerX = ceil(Nx / 2);
centerY = ceil(Ny / 2);
theta = atan2(Y - centerY, X - centerX); % 角度
rho = sqrt((X - centerX).^2 + (Y - centerY).^2); % 半径

% 定义角度分辨率
numBins = 360; % 每度一个bin
thetaBins = linspace(-pi, pi, numBins+1);

% 初始化存储每个角度的能量
angleEnergy = zeros(numBins, 1);

% 计算每个角度范围的能量
for i = 1:numBins
    mask = theta >= thetaBins(i) & theta < thetaBins(i+1);
    angleEnergy(i) = sum(energyDistribution(mask));
end

% 归一化为百分比
angleEnergyPercentage = (angleEnergy / totalEnergy) * 100;

% 绘制角度能量分布
figure;
polarplot(thetaBins(1:end-1), angleEnergyPercentage, 'r-');
title('Angular Energy Distribution (%)');
ax = gca;
ax.RLim = ;

### 说明：

1. **极坐标转换**：
   - 使用 `atan2` 计算每个点的角度 `theta`，以中心 `(centerX, centerY)` 为基点。
   - 半径 `rho` 用于对能量进行角度与半径的区分。

2. **角度能量计算**：
   - 将空间频率数据按角度划分为多个段。
   - 使用掩码 `mask` 提取对应角度段内的能量值并求和。

3. **结果可视化**：
   - `polarplot` 绘制不同角度下的能量分布百分比，直观展示了频域的角度分布特性。

通过以上步骤，你可以查看不同角度所占据的总能量比例，并发现频域中可能存在的非对称性或集中特征。如果有其他问题或需求，欢迎随时提出！
--- 光学专家Gpt

Anonymous 发表于 2025-1-6 10:37:45

我觉得这个角度分布不太对，或者你直接作图关于频数和能量的分布吧，展现频域的能量分布特性

机器人Gpt 发表于 2025-1-6 10:37:57

LIYU 发表于 2025-01-06 10:37
我觉得这个角度分布不太对，或者你直接作图关于频数和 ...

为了展示频域中的能量分布特性，我们可以通过对傅里叶变换结果的频数进行径向平均，并绘制其能量分布。这样可以直观地看到不同频率下的能量占比。

以下代码展示了如何实现这一点：

% 打开 ZBF 文件
filename = 'your_file.zbf'; % 请替换为实际的文件名
fileID = fopen(filename, 'r', 'l'); % 'l' denotes little-endian

if fileID == -1
    error('无法打开文件');
end

% 假设我们知道数据大小为 Nx x Ny
Nx = 256; % 从文件中解析
Ny = 256; % 从文件中解析
dataSize = Nx * Ny;

% 读取照度数据
intensity = fread(fileID, dataSize, 'double');

% 读取相位数据
phase = fread(fileID, dataSize, 'double');

fclose(fileID);

% 对照度和相位数据进行预处理
intensity(isnan(intensity) | isinf(intensity)) = 0;
phase(isnan(phase) | isinf(phase)) = 0;

% 对照度数据进行归一化处理
intensity = intensity / max(intensity(:));

% 确保数据组织成矩阵形式
intensity = reshape(intensity, Nx, Ny);
phase = reshape(phase, Nx, Ny);

% 形成复数形式的光场
complexField = sqrt(intensity) .* exp(1i * phase);

% 修正可能的 NaNs/Infs 后组合
complexField(isnan(complexField) | isinf(complexField)) = 0;

% 傅里叶变换
fftResult = fftshift(fft2(complexField));

% 能量分布
energyDistribution = abs(fftResult).^2;

% 创建网格来计算半径
= meshgrid(1:Nx, 1:Ny);
centerX = ceil(Nx / 2);
centerY = ceil(Ny / 2);
rho = sqrt((X - centerX).^2 + (Y - centerY).^2); % 半径

% 最大半径
maxRadius = floor(sqrt(centerX^2 + centerY^2));

% 初始化存储每个半径的能量
radialEnergy = zeros(maxRadius, 1);
count = zeros(maxRadius, 1);

% 计算径向能量分布
for x = 1:Nx
    for y = 1:Ny
      radius = round(rho(y, x));
      if radius <= maxRadius
            radialEnergy(radius) = radialEnergy(radius) + energyDistribution(y, x);
            count(radius) = count(radius) + 1;
      end
    end
end

% 计算平均能量
radialEnergy(count > 0) = radialEnergy(count > 0) ./ count(count > 0);

% 归一化为百分比
totalEnergy = sum(radialEnergy);
radialEnergyPercentage = (radialEnergy / totalEnergy) * 100;

% 绘制径向能量分布
figure;
plot(1:maxRadius, radialEnergyPercentage, 'b-');
xlabel('Frequency Radius');
ylabel('Energy Distribution (%)');
title('Radial Energy Distribution in Frequency Domain');

### 说明：
1. **径向平均**：通过计算距离中心的半径 `rho`，对每个半径距离内的能量进行累积并求平均。
2. **归一化**：计算每个径向距离上的能量占总能量的百分比，以便更清楚地表示频域的能量分布特性。
3. **结果可视化**：使用 `plot` 函数绘制频数（以半径为单位）与能量分布的关系图。

这份代码将帮助你更好地理解原始信号在不同频率成分上的能量分布特性。根据实际需求，你可以进一步调整代码的各个部分。
--- 光学专家Gpt

Anonymous 发表于 2025-1-6 10:42:10

数组索引必须为正整数或逻辑值。

出错 read_ZBF (line 63)
            radialEnergy(radius) = radialEnergy(radius) + energyDistribution(y, x);

机器人Gpt 发表于 2025-1-6 10:42:24

LIYU 发表于 2025-01-06 10:42
数组索引必须为正整数或逻辑值。

出错 read_ZBF (lin ...

这个错误提示是因为在试图索引数组时，使用了一个非正整数值。在这个上下文中，`radius` 的计算结果应该是一个正整数，但可能存在计算误差或未考虑的边缘情况，使得 `radius` 变为零或其他不能用作索引的值。

为了解决这个问题，我们需要确保 `radius` 在索引数组时是一个有效的正整数。通常，我们会将 `radius` 的计算结果限定在一个合理的范围内，这样可以防止越界或无效索引。

以下是修改后的代码：

% 打开 ZBF 文件
filename = 'your_file.zbf'; % 请替换为实际的文件名
fileID = fopen(filename, 'r', 'l'); % 'l' denotes little-endian

if fileID == -1
    error('无法打开文件');
end

% 假设我们知道数据大小为 Nx x Ny
Nx = 256; % 从文件中解析
Ny = 256; % 从文件中解析
dataSize = Nx * Ny;

% 读取照度数据
intensity = fread(fileID, dataSize, 'double');

% 读取相位数据
phase = fread(fileID, dataSize, 'double');

fclose(fileID);

% 对照度和相位数据进行预处理
intensity(isnan(intensity) | isinf(intensity)) = 0;
phase(isnan(phase) | isinf(phase)) = 0;

% 对照度数据进行归一化处理
intensity = intensity / max(intensity(:));

% 确保数据组织成矩阵形式
intensity = reshape(intensity, Ny, Nx); % 确保维度顺序正确
phase = reshape(phase, Ny, Nx);

% 形成复数形式的光场
complexField = sqrt(intensity) .* exp(1i * phase);

% 修正可能的 NaNs/Infs 后组合
complexField(isnan(complexField) | isinf(complexField)) = 0;

% 傅里叶变换
fftResult = fftshift(fft2(complexField));

% 能量分布
energyDistribution = abs(fftResult).^2;

% 创建网格来计算半径
= meshgrid(1:Nx, 1:Ny);
centerX = ceil(Nx / 2);
centerY = ceil(Ny / 2);
rho = sqrt((X - centerX).^2 + (Y - centerY).^2); % 半径

% 最大半径
maxRadius = floor(sqrt(centerX^2 + centerY^2));

% 初始化存储每个半径的能量
radialEnergy = zeros(maxRadius, 1);
count = zeros(maxRadius, 1);

% 计算径向能量分布
for x = 1:Nx
    for y = 1:Ny
      radius = round(rho(y, x));
      if radius >= 1 && radius <= maxRadius
            radialEnergy(radius) = radialEnergy(radius) + energyDistribution(y, x);
            count(radius) = count(radius) + 1;
      end
    end
end

% 平均化每个半径的能量
radialEnergy(count > 0) = radialEnergy(count > 0) ./ count(count > 0);

% 归一化为百分比
totalEnergy = sum(radialEnergy);
radialEnergyPercentage = (radialEnergy / totalEnergy) * 100;

% 绘制径向能量分布
figure;
plot(1:maxRadius, radialEnergyPercentage, 'b-');
xlabel('Frequency Radius');
ylabel('Energy Distribution (%)');
title('Radial Energy Distribution in Frequency Domain');

### 说明：

1. **维度确保**：在 reshaping intensity 和 phase 数据时，确保它们的顺序对齐 `Ny` 和 `Nx`。
2. **有效索引**：添加条件检查，确保 `radius` 在1到`maxRadius`之间。
3. **网格生成**：修正生成网格的代码，确保 X 和 Y 被定义，以便用于半径计算。

这些调整应该解决之前遇到的问题，并为你提供所需的频域能量分布特性。
--- 光学专家Gpt

查看完整版本: ZBF