通信算法之169：通信物理层QPSK中的FFT估计载波频偏算法

一. 求频偏几种常见的算法

1. 训练序列相关

2. FFT算法

3. 锁相环

二. FFT算法

       快速傅里叶变换（FFT）的频偏估计算法，它能直接检测得到载波频偏绝对值大小，且精度与输入端信噪比无关，可应用于中继卫星通信中高速数传且频带有效的同步解调。

三.库函数代码

%% Compensate for Frequency Offset in a QPSK Signal
% Compensate for a 4 kHz frequency offset imposed on a noisy QPSK signal.
% 
% Set the example parameters.

nSym = 2048;       % Number of input symbols
sps = 4;           % Samples per symbol
nSamp = nSym*sps;  % Number of samples
fs = 80000;        % Sampling frequency (Hz)
%% 
% Create a square root raised cosine transmit filter.
%%
txfilter = comm.RaisedCosineTransmitFilter(...
    'RolloffFactor',0.2, ...
    'FilterSpanInSymbols',8,   ...
    'OutputSamplesPerSymbol',sps);
%% 
% Create a phase frequency offset object to introduce the 4 kHz frequency 
% offset.
%%
freqOffset = comm.PhaseFrequencyOffset(...
    'FrequencyOffset',-4000, ...
    'SampleRate',fs);
%% 
% Create a coarse frequency compensator object to compensate for the offset.
%%
freqComp = comm.CoarseFrequencyCompensator(...
    'Modulation','QPSK', ...
    'SampleRate',fs, ...
    'FrequencyResolution',1);
%% 
% Generate QPSK symbols, filter the modulated data, pass the signal through 
% an AWGN channel, and apply the frequency offset.
%%
data = randi([0 3],nSym,1);
modData = pskmod(data,4,pi/4);
txSig = txfilter(modData);
rxSig = awgn(txSig,20,'measured');
offsetData = freqOffset(rxSig);
%% 
% Compensate for the frequency offset using |freqComp|. When the frequency 
% offset is high, it is beneficial to do coarse frequency compensation prior to 
% receive filtering because filtering suppresses energy in the useful spectrum.
%%
[compensatedData,estFreqOffset] = freqComp(offsetData);
%% 
% Display the estimate of the frequency offset.
%%
estFreqOffset
%% 
% Return information about the |freqComp| object. To obtain the FFT length, 
% you must call |freqComp| prior to calling the |info| method.
%%
freqCompInfo = info(freqComp)
%% 
% Create a spectrum analyzer object and plot the offset and compensated 
% spectra. Verify that the compensated signal has a center frequency at 0 Hz and 
% that the offset signal has a center frequency at -4 kHz.
%%
specAnal = dsp.SpectrumAnalyzer('SampleRate',fs,'ShowLegend',true, ...
    'ChannelNames',{'Offset Signal' 'Compensated Signal'});
specAnal([offsetData compensatedData])

四.工程实现代码

% QPSK FFT方法
% QQ : 1279682290

clc
clear
close all

%-参数---------------------------
sps = 4;%-升余弦滚降滤波器
Fs = 1000*sps; %采样频率

% Fs = 1000; %采样频率
T = 1/Fs;  %采样周期
L = 1500;
t = (0:L-1)*T;
offset = 100;

data = randi([0 3],L,1);
qpsk_symbol0 = pskmod(data,4,pi/4).';
% qpsk_symbol0 = exp(1j*(randi([0 3],1,L)*pi/2 + pi/4));

scatterplot(qpsk_symbol0)

%-升余弦滚降滤波器-------------------------------
txfilter = RaisedCosineTnsmitFilter(...
    'RolloffFactor',0.25, ...
    'FilterSpanInSymbols',8,   ...
    'OutputSamplesPerSymbol',sps);
qpsk_symbol = txfilter(qpsk_symbol0.').';

%-----------------------------------------
% spectrumAnalyzer = dsp.SpectrumAnalyzer('SampleRate',Fs);
% spectrumAnalyzer(qpsk_symbol)

%--ADD 频偏---------------------------------------
freqOffset = PhaseFreqncyOffset(...
    'FrequencyOffset',offset, ...
    'SampleRate',Fs);
S = freqOffset(qpsk_symbol.').';


%--ADD 频偏---------------------------------------

S_addnoise = awgn(S,20,'measured');

L = length(qpsk_symbol);


ylabel('幅度');

figure;
Y2 = fft( );
x_fre = (0:L-1)./L*Fs;
plot(x_fre,abs(Y2/L));
xlabel('频率');
ylabel('幅度');
title('无噪声时的信号频率图')

% 方法一：
figure;
Y3 = fft( );
x_fre = (0:L-1)./L*Fs;
plot(x_fre,abs(Y3/L));
xlabel('频率');
ylabel('幅度');
title('有噪声时的信号频率图')

[maxvaule,index] = max(Y3);
foffset_esti = x_fre(1,index)/4;

%方法二：FFT方法
freqComp = arseFreqcyCompensator(...
    'Modulation','QPSK', ...
    'SampleRate',Fs, ...
    'FrequencyResolution',1);
[compensatedData,estFreqOffset] = freqComp(S_addnoise.');

FFT估计频偏的方法使得频偏偏差已经比较小时，随着时间的积累依然会造成星座图相位旋转，这样会对解调造成误判，尝试用以下方法解决：

（1）频偏比较小后，进入锁相环；这个我之前已经仿真过BPSK和QPSK采用锁相环来纠正固定频偏和相偏的误码率曲线，和理论曲线相差不大。当然，调锁相环的参数方面经验尚且不足。

（2）如果信号本身已经采用差分编码，对频偏也能有一定程度的容忍。比如采用DBPSK或者DQPSK，因为相邻符号的相差用来做判决。

（3）当帧长持续时间Tframe较短，在Tframe时间内，频偏偏差不会造成QPSK星座图旋转超过pi/4，也不会对判决有影响，等等。

五.升余弦滤波器

作用

升余弦滤波器常常作为发送端的成型滤波器，用来抑制信号带外辐射。

MATLAB实现

在maltab中，主要用rcosdesign函数来实现FIR脉冲成型滤波器，其他相关的函数还有comm.RaisedCosineTransmitFilter，comm.RaisedCosineTransmitFilter。

语法：

b = rcosdesign(beta,span,sps)

b = rcosdesign(beta,span,sps,shape)

描述：

b = rcosdesign(beta,span,sps)返回系数b，对应于一个滚降系数为beta的均方根升余弦FIR滤波器。滤波器截断为span个symbols，每个symbol包含sps个采样点。滤波器的阶数sps*span必须是偶数。滤波器能量为1。

b = rcosdesign(beta,span,sps，shape)返回一个均方根升余弦滤波器（shape为’sqrt’）或者一个升余弦滤波器（shape为’normal’）。

升余弦滤波器的例子

升余弦可以抑制码间干扰，发送和接收滤波器分开。

升余弦滤波器的主要参数是其滚降系数，滚降系数直接决定了滤波器的带宽。理想的升余弦滤波器由无限多个抽头。因此，实际上的升余弦滤波器是加窗的。窗的长度由FilterSpanInSymbols特性来控制。在本例中，我们将窗长度设定为6个symbol长度。这样的滤波器也有3个symbol的群延时。升余弦滤波器用在信号上采样后的脉冲成型。因此，我们也需要指定上采样倍数。升余弦滤波器的参数如下：

Nsym=6;

beta=0.5;

samsPerSym=8;

【参考：MATLAB分析升余弦滤波器_sps samples symbol_Brandon懂你的博客-CSDN博客】

【参考：知网论文《快速傅里叶变换载波频偏估计算法》】

用升余弦发射滤波器系统结构体构建滤波器并用fvtool来可视化滤波器特性。

这个结构体设计了一个直接型多项式单位能量FIR滤波器。滤波器阶数Nsym*sampsPerSym，抽头数Nsym*sampsPerSym+1。可以利用gain特性归一化滤波器系数以使滤波后的和未滤波数据匹配。