利用ECEF坐标和仰角计算多普勒频移

要计算多普勒频移，可以使用以下公式：

$$f_D = -\frac{2v}{c}f_0\sin\theta$$

其中，$v$ 是接收器和发射器之间的相对速度，$c$ 是光速，$f_0$ 是发射频率，$\theta$ 是信号的仰角。

如果已知接收器和发射器的 ECEF 坐标，可以使用以下方法计算相对速度 $v$ 和信号的仰角 $\theta$：

% 输入参数
f0 = 10e9;                  % 发射频率
rE = 6378137;               % 地球半径
c = 3e8;                    % 光速
receiverECEF = [0, 0, 0];   % 接收器 ECEF 坐标，单位为米
satelliteECEF = [1e7, 1e7, 1e7]; % 卫星 ECEF 坐标，单位为米
receiverLLH = ecef2lla(receiverECEF);   % 将接收器 ECEF 坐标转换为经纬高
satelliteLLH = ecef2lla(satelliteECEF); % 将卫星 ECEF 坐标转换为经纬高

% 计算相对速度
receiverENU = ecef2enu(receiverECEF, receiverLLH);    % 将接收器 ECEF 坐标转换为东北天坐标系
satelliteENU = ecef2enu(satelliteECEF, receiverLLH);  % 将卫星 ECEF 坐标转换为接收器所在经度的东北天坐标系
vENU = satelliteENU - receiverENU;                   % 计算卫星和接收器之间的相对速度（东北天坐标系）
vECEF = enu2ecef(vENU, receiverLLH);                 % 将相对速度转换为 ECEF 坐标系
v = norm(vECEF);                                    % 计算相对速度的模长

% 计算仰角
[az, el, ~] = geodetic2aer(satelliteLLH(1), satelliteLLH(2), satelliteLLH(3), receiverLLH(1), receiverLLH(2), receiverLLH(3), wgs84Ellipsoid);
theta = pi/2 - el;  % 将仰角转换为弧度制，并计算相对于垂直方向的角度

% 计算多普勒频移
fD = -(2*v/c)*f0*sin(theta);

% 输出结果
fprintf('多普勒频移为: %.2f Hz\n', fD);

请注意，此代码中使用了 MATLAB 自带的一些函数来进行 ECEF 坐标系和经纬高坐标系之间的转换，包括 `ecef2lla`、`ecef2enu`、`enu2ecef`、`geodetic2aer`。