代码片记录-使用DDC对两路信号相位差进行求取（matlab实现）

参考链接：

1、FPGA综合系统设计（七）基于DDC的两路信号相位差检测_fpga相位测量-CSDN博客 （真的太牛了）

2、双通道中频信号数字下变频及相位差估计（FPGA）_根据 i、q 两路信号以如何求出中频信号的频率和相位差-CSDN博客 （通过这个篇博客注意到了反正切的值域问题，里面其实还有fapa实现的具体实现，之后可以参考）

3、atan2函数的原理 - 百度文库 (baidu.com)  (反正切函数的相位矫正)

4、相位差_百度百科 (baidu.com)                   （相位差的定义，常识）

matlab的实现在F:\WorkSpace\practice\1、D9射频电源项目20240118-\我自己的之前找的资料\ADC采样\波形复原\DDC_WITH_PHASE_ABTAIN_REMOUDLE_DEMO_V1_0.m

思考：

1、使用iq解调的时候，iq的选取，不会影响幅值的求取，但是对相位的求取会有一定的影响，下面进行讨论：

• A 、使用直接iq解调，在一路信号上求取幅值
• 

• 最后求取的幅值都是Aif，但是这里求得的相位是中频的初始相位，是不会变化的，他和原始信号有个固定的关系

• B、使用ddc的解调方式，选取两路作为iq信息
• 

• 求得的也是幅值，但是相位是一个和时间有关的函数，是中频信号的实际相位，和时间有关的

最开始我的疑问点：

1、使用一路DDC求取所求取得到的相位信息是什么？

是混频之后的中频的真实（实时）相位，这个和时间有关，是一个中频频率和时间的函数，初始相位是原始信号和本征信号的相位差。

2、如何获得两路（电流电压）信号的相位差，将两个信号都进行下变频，将最后的相位信息相减，得到的就是两路信号的相位差，下面是推导过程

3、采用matlab进行仿真模拟

给出一个模拟的波形，电流和电压的相位相差40度，最终运算结果也是40度，有了这个数据之后可以计算驻波比，以及其它的有关阻抗匹配的参数了

原始信号存在相位差

通过计算得到相位差结果，比较稳定，原理参照上面的手写笔记

4、出现了一个问题，我好像都不知道着相位的定义了，里面并没有夹杂rf的实时相位呀，为啥会不到90，有一个一直是80的样子，有点奇怪哇，

问题找到了是反正切函数的值域问题，用的是两象限的反正切函数，值域是-0.5pi到0.5pi，需要对其进行矫正，但是我尝试了手动按照象限去矫正，相位的波形畸变严重，有没有其它的方法？

找到了修正方法，真的蠢好吧，之后还是可以按照这个实现的

4、程序记录文章来源地址https://www.toymoban.com/news/detail-850554.html

%% 描述
% V1.0--------基于IQ解调的说明
%       用于测试IQ正交解调，求幅值信息
%       使用了fir（汉宁窗口）滤波器对混频信号进行高频滤除
%       后记采用了滑动平均滤波器进行幅值输出
% V2.0------- 基于原始信号相位差求取实现
%       对两路原始信号分别进行DDC处理，得到相位信息之后，将两路相位相减，得到原始信号相位差啊
%       仿真验证这个是可行的，需要注意的点是反正切函数的值域问题，最开始就是因为这个没注意，然后走了很多弯路
%% 开始
close all;         %先关闭所有图片
clear all;

SW=1;              %信号源控制开关             0--lia 抓取数据   1--理想数据
SW1=1;             %数字滤波方法选择开关       0--平均值滤波     1--滑动窗口滤波 
SW2=1;             %低通滤波器选择开关         0--未知低通滤波   1--fir滤波器
SW3=0;             %数据突变控制开关           0--不开启数据突变 1--开启数据突变
%% 采样数据初始化
Fs=125000000;      %采样频率(Hz)
Ft=40680000;       %模拟信号频率 
N=511;             %采样点数
real_nums=N+1;
t=[0:1/Fs:N/Fs];   %采样时刻,表示采了多少个数据点
Fm=30000000;       %本征频率 30MHz
A1=1;              %本征信号的幅值
P1=0;              %本征信号相位 角度
P3=40;             %原始信号相位差定义 角度制
count=10;          %用于均值处理
for_nums=10;       %for循环次数
amp_filter=eye(1,real_nums); %存储滤波之后幅值数据

%% 本征信号生成 
sm1=A1*cos(2*pi*Fm*t+pi*P1/180);
sm2=-A1*sin(2*pi*Fm*t+pi*P1/180);

%% adc数据获取
if SW==0
sin_wave=eye(1,real_nums);
sheetname='sheet1';
filename='LIA采集到的adc数据用于matlab仿真fft.csv';
% Data=readtable(filename,'Sheet',sheetname);
data=xlsread(filename,'LIA采集到的adc数据用于matlab仿真fft','F2:F513');
sine_wave1(1:real_nums)=data';
end
if SW==1
    sine_wave1=600*cos(2*pi*Ft*t);
    sine_wave2=600*cos(2*pi*Ft*t+pi*P3/180);
end

%% 采样值突变
if SW3==1
    for i=200:210
        sine_wave1(i)=sine_wave1(i)*5*(i-200)/10;
    end
end

%% 显示原始信号
figure(1);
subplot(3,1,1);
plot(t,sine_wave1);
title('“电压”原始信号');
subplot(3,1,2);
plot(t,sine_wave2);
title('”电流“原始信号');
subplot(3,1,3);
plot(t,sm1);
title('本征信号');

%% IQ 正交处理
figure(2);
wave_of_deal_i1=sm1.*sine_wave1;
subplot(4,1,1);
plot(t,wave_of_deal_i1);
title('I1信号');

wave_of_deal_q1=sm2.*sine_wave1;
subplot(4,1,2);
plot(t,wave_of_deal_q1);
title('Q1信号');

wave_of_deal_i2=sm1.*sine_wave2;
subplot(4,1,3);
plot(t,wave_of_deal_i2);
title('I2信号');

wave_of_deal_q2=sm2.*sine_wave2;
subplot(4,1,4);
plot(t,wave_of_deal_q2);
title('I2信号');

%% 低通滤波
%   ----------------- 未知低通滤波方法
if SW2==0
result_wave_i=lowpass(wave_of_deal_i,10e6,30e6);
result_wave_q=lowpass(wave_of_deal_q,10e6,30e6);
end
%   ----------------- fir 低通滤波
if SW2==1
FilterLen = 512;          % 滤波器长度
FilterStep = FilterLen-1; % 滤波器阶数
FreqEnd  =15e6;           % 高截止频率
w2 = FreqEnd /Fs*2*pi;    % 归一化
window=hann(FilterLen);   % 汉宁窗，使得fir滤波器第一个旁瓣增益降低，改善滤波器性能
%window=rectwin(FilterLen); % 矩形窗，窗口区间都是为1

hn=fir1(FilterStep,w2/(pi),'low',window);                     % 生成低通滤波器，暂时没有使用，被滤波器设计器的滤波器取代
figure(3);
freqz(hn,1,512);
hm=hann_filter;
result_wave_i1 = conv(wave_of_deal_i1,hm.Numerator);          % 对原始序列和滤波器序列求卷积得到滤波后的结果
result_wave_q1 = conv(wave_of_deal_q1,hm.Numerator);          % 这里直接使用了滤波器设计工具生成的滤波器，但是都是hn窗的
result_wave_i1 = result_wave_i1(FilterLen/2:FilterLen/2+N);   % 滤波结果是居中的数据
result_wave_q1 = result_wave_q1(FilterLen/2:FilterLen/2+N);   % 滤波结果是居中的数据

result_wave_i2 = conv(wave_of_deal_i2,hm.Numerator);          % 对原始序列和滤波器序列求卷积得到滤波后的结果
result_wave_q2 = conv(wave_of_deal_q2,hm.Numerator);          % 这里直接使用了滤波器设计工具生成的滤波器，但是都是hn窗的
result_wave_i2 = result_wave_i2(FilterLen/2:FilterLen/2+N);   % 滤波结果是居中的数据
result_wave_q2 = result_wave_q2(FilterLen/2:FilterLen/2+N);   % 滤波结果是居中的数据
end

figure(4);
subplot(4,1,1);
plot(t,result_wave_i1);
title('i1低通处理');

subplot(4,1,2);
plot(t,result_wave_q1);
title('q1低通处理');

subplot(4,1,3);
plot(t,result_wave_i2);
title('i2低通处理');

subplot(4,1,4);
plot(t,result_wave_q2);
title('q2低通处理');
%% 幅值求取
amp=sqrt(result_wave_i1.*result_wave_i1+result_wave_q1.*result_wave_q1);
figure(5);
subplot(2,1,1);
plot(t,amp*2);   
title('幅值结果');

%% 均值滤波        每20个点求取一次均值
if SW1==0
    while count<real_nums
        amp_temp=0;
        for k=1:for_nums
            amp_temp=amp_temp+amp(count+k-for_nums);
        end
        count=count+1;
        amp_temp=amp_temp/for_nums;
        amp_filter(count)=amp_temp;
    end
end

%% 滑动窗口滤波
if SW1==1
    windowsize=10;
    b=(1/windowsize)*ones(1,windowsize);
    a=1;
    amp_filter=filter(b,a,amp);
end

subplot(2,1,2);
plot(t,amp_filter*2);   
title('滤波之后幅值');

%% 相位求取
figure(6);
% phase1=atan2(result_wave_q1,result_wave_i1)*180/pi;      %求四象限上的反正切，取值范围是-pi到pi
%---------------求两个象限上的反正切， 手动修正，效果相同--------------------%
phase1=atan(result_wave_q1./result_wave_i1)*180/pi;
%phase1的修正
for index=1:1:real_nums
    if and((result_wave_i1(index)>=0),(result_wave_q1(index)>=0)) %第一象限
       phase1(index)=phase1(index);
    end
    if and(result_wave_i1(index)<0,result_wave_q1(index)>=0)     %第二象限
       phase1(index)=phase1(index)+180;
    end
    if and(result_wave_i1(index)<0,result_wave_q1(index)<0 )     %第三象限
       phase1(index)=phase1(index)-180;
    end
    if and(result_wave_i1(index)>=0,result_wave_q1(index)<0)     %第四象限
       phase1(index)=phase1(index);
    end
end
%----------------------------end------------------------------------------%
subplot(3,1,1);
plot(t,phase1);
title('电压中频相位');

%  phase2=atan2(result_wave_q2,result_wave_i2)*180/pi;     %求四象限上的反正切，取值范围是-pi到pi
 %---------------求两个象限上的反正切， 手动修正，效果相同-------------------%
phase2=atan(result_wave_q2./result_wave_i2)*180/pi;
%phase2的修正
for index=1:1:real_nums
    if and((result_wave_i2(index)>=0),(result_wave_q2(index)>=0))%第一象限
       phase2(index)=abs(phase2(index));
    end
    if and(result_wave_i2(index)<0,result_wave_q2(index)>=0)     %第二象限
       phase2(index)=phase2(index)+180;
    end
    if and(result_wave_i2(index)<0,result_wave_q2(index)<0 )     %第三象限
       phase2(index)=phase2(index)-180;
    end
    if and(result_wave_i2(index)>=0,result_wave_q2(index)<0)     %第四象限
       phase2(index)=phase2(index);
    end
end
%----------------------------end------------------------------------------%
subplot(3,1,2);
plot(t,phase2);
title('电流中频相位');

detal=phase2-phase1;

% 直接求相位会出现错误，需要对计算相位进行修正
% 修正原则是当啊相位差超过pi时减去2pi，小于-pi时，加上2pi
for index=1:1:real_nums
    if detal(index)>180
        detal(index)= detal(index)-360;
    end
    if detal(index)<-180
        detal(index)= detal(index)+360; 
    end
end

subplot(3,1,3);
plot(t,detal);
title('原始信号相位差');


%% end