作者：桂。
时间：2017-12-02 23:29:48
链接：http://www.cnblogs.com/xingshansi/p/7956491.html
一、相位提取
以正弦信号为例，x = sin(2pi*f*t+pi)，希望提取phi：
思路1：通过Hilbert变化解决
思路2：借助FFT，利用插值思想，估计Phi;
思路3：借助全相FFT（apFFT, all phase FFT）实现。
思路三可提取信号相位，这一点FFT做不到，而相位信息通常可判断相位调制类型，可用于情报的脉内检测。
全相FFT思路：
- 选定N点窗，如hanning
- 窗函数自相关，并归一化
- 对2N-1序列x(n)加窗，
- 将2N-1个点，每间隔N点相加
- FFT实现
二、仿真验证
clc;clear all;close all;fs = 1e9; fo = 200e6; t = 0:1/fs:1023/fs; tao = 0.3/3e8*sind(45); SNR = 20; ch1 = awgn(sin(2*pi*t*fo),SNR) ; ch2 = awgn(sin(2*pi*t*fo + 2*pi*tao*fo),SNR); % ch1 = sin(2*pi*t*fo) ; % ch2 = sin(2*pi*t*fo + 0.5); pha = angle(hilbert(ch2))-angle(hilbert(ch1)); figure() subplot 211 plot(t*fs,pha); subplot 212 plot(t(1:512)*fs,abs(fft(ch1(1:512))),'r--');hold on; %FFT提取相位 pha1 = angle(fft(ch1(1:512)).*fft(ch2(1:512))); figure() subplot 211 plot(t(1:512)*fs,pha1); subplot 212 plot(t(1:512)*fs,abs(fft(ch1(1:512))),'r--');hold on; %apFFT提取相位 win = hanning(512)'; win1 = conv(win,win); win1 = win1/sum(win1); y1 = ch1(1:1023).*win1; y2 = ch2(1:1023).*win1; out1 = [0,y1(1:511)]+y1(512:1023); out2 = [0,y2(1:511)]+y2(512:1023); pha2 = angle(fft(out1).*conj(fft(out2))); figure() subplot 211 plot(t(1:512)*fs,pha2); subplot 212 plot(t(1:512)*fs,abs(fft(ch1(1:512))),'r--');hold on; [~,pos] = max(abs(fft(ch1(1:512)))); [pha2(pos) mean(pha) ;-pi+ 2*pi*tao*fo 2*pi*tao*fo] theta_est = asind((pha2(pos))/2/pi/fo/0.3*3e8)+90; abs(theta_est-45)
另外，频谱细化，可借助zoom-FFT：
fs = 2048; T = 100; t = 0:1/fs:T; x = 30 * cos(2*pi*110.*t) + 30 * cos(2*pi*111.45.*t) + 25*cos(2*pi*112.3*t) + 48*cos(2*pi*113.8.*t)+50*cos(2*pi*114.5.*t); [f, y] = zfft(x, 109, 115, fs); plot(f, y);function [f, y] = zfft(x, fi, fa, fs) % x为采集的数据 % fi为分析的起始频率 % fa为分析的截止频率 % fs为采集数据的采样频率 % f为输出的频率序列 % y为输出的幅值序列(实数）f0 = (fi + fa) / 2; %中心频率 N = length(x); %数据长度r = 0:N-1; b = 2*pi*f0.*r ./ fs; x1 = x .* exp(-1j .* b); %移频bw = fa - fi; B = fir1(32, bw / fs); %滤波 截止频率为0.5bw x2 = filter(B, 1, x1); c = x2(1:floor(fs/bw):N); %重新采样 N1 = length(c); f = linspace(fi, fa, N1); y = abs(fft(c)) ./ N1 * 2; y = circshift(y, [0, floor(N1/2)]); %将负半轴的幅值移过来 end
参考：Digital Receiver-based Electronic Intelligence System Configuration for the Detection and Identification of Intrapulse Modulated Radar Signals
)
