激光多普勒测振仿真，反正切法与微分交叉相乘算法

本文仿真了一个微振动引起的多普勒频移，通过多普勒频移还原物体振动信号

原理为：多普勒频移fd = 4*pi*st/激光波长

通过将信号分为IQ两路后，通过反正切或微分交叉相乘就能得到初始信号

% 研究噪声对信号解调的影响
% 这段代码未加白噪声
%用awgn添加高斯白比例为34、22、12、2db

close all;
clear;
clc;

fs=400e6;%采样率
dt=1/fs;%采样间隔
T=2e-2;%采样时间
t=0:dt:T-dt;
[M,N]=size(t);
Fr=40e6;%中频
laser = 1550e-9; %激光波长

fv=4000;%简谐振动频率
A=1e-7;%简谐振动振幅
Sv=A.*cos(2*pi*fv*t);%简谐振动信号
% Vv=-A.*2*pi*fv*sin(2*pi*fv*t);%简谐振动速度
% Fv=2*Vv/1.55e-6;
% max(Vv)
% min(Vv)
% max(Fv)
% min(Fv)


S=cos(2*pi*Fr*t+4*pi*Sv/laser);%振动表现为相位的改变
% S=awgn(S,10);%添加高斯白噪声
% [spec,f] = myFFT(S,fs);
% plot(f,spec)
figure;
plot(t,S);

S1 = filter(filter_high3,S);
B = fft(S1);
B = B(1:N/2+1);
f = (0:N/2)*fs/N;
figure(1);
plot(t,S1);title('滤波后的采样信号')
figure(2);
plot(f,abs(B));title('频谱图')
xlim([3.999e7 4.001e7])
figure;
S_db=20*log10(abs(B));
plot(f,S_db);
xlim([3.999e7 4.001e7])

xlabel('频率(Hz)');
ylabel('幅值(db)');


L1 = length(S);
t1 = (0:L1-1)/fs;
I=cos(2*pi*Fr*t).*S;
Q=sin(2*pi*Fr*t).*S;


%低通滤波器设计
%滤波后的数据
It = filter(lowfilter1,I);
Qt = filter(lowfilter1,Q);
figure;
subplot(221),plot(t1,I),title('I路信号');
xlim([0.5e-3 2e-3]);
subplot(222),plot(t1,It),title('滤波后I路信号');
subplot(223),plot(t1,Q),title('Q路信号');
xlim([0.5e-3 2e-3]);
subplot(224),plot(t1,Qt),title('滤波后Q路信号');


%方法一：相位反正切
p =atan(Qt./It);
p1 = p;
%相位解包裹
for i = 2:length(p)
    diff = p(i) - p(i-1);
    if diff > pi/2
        p1(i:end) = p1(i:end) - pi;
    elseif diff < -pi/2
        p1(i:end) = p1(i:end) + pi;
    end
end

% 微分交叉相乘算法
% x11 = diff(It);
% x22 = diff(Qt);
% L = length(It);
% x_d = (It(1:L-1)).*x22-(Qt(1:L-1)).*x11;

% X = zeros(1,L-1);
% 积分
% for j = 1:L-1
%    X(j) = 0;
% end
% for i = 2:L-1
%     X(i) = X(i-1)+(x_d(i))/(It(i).^2+Qt(i).^2);
% end

st = laser/(4*pi).*p1;
figure;
plot(st);
hold on;
Sv1=A.*cos(2*pi*fv*t+pi/2);
plot(Sv1);
legend('解调信号','实际信号');