\(\begin{align*} \newcommand{\dif}{\mathop{}\!\mathrm{d}} \newcommand{\belowarrow}[1]{\mathop{#1}\limits_{\uparrow}} \newcommand{\bd}{\boldsymbol} \newcommand{\L}{\mathscr{L}} \newcommand{\xleftrightarrow}[1]{\stackrel{#1}{\longleftrightarrow}} \end{align*}\)
实验要求
- 验证性实验
- 用 Matlab 验证DFT运算的对称性质
- 编程计算出一个2N点的实序列的频谱,要求只使用一次 N 点的 FFT 和其他运算,并与一次 2N 点 FFT 作出的频谱进行比较验证
- 应用性实验
- 去中国期刊网、IEEE网站等查找频率估计方法,利用至少两种频率估计方法编程求给定信号的频率 ,要求:
- 所有算法使用相同点数的FFT(小于等于100点)
- 无噪声、20dB、10dB、5dB、0dB、-5dB,每个信噪比至少100次(无噪声除外)
- 绘出估计均方误差图
- 信号来源:用采样率8000Hz,对频率为350Hz的单频正弦信号进行采样,取70个点作为原始信号(无噪声)
- 线性调频信号LFM: $x(nT)=\exp(j2\pi f_0 nT+j \pi u (nT)^2)$,$f_0$ 为中心频率,$T$ 为信号采样间隔,$u$ 为调频率,要求:
- 在LFM信号中加入噪声信号 $n(t)$,采用匹配滤波器的方法对LFM信号进行脉冲压缩,并观察效果。$h(t)=s^*(t_0-t)$
- 查阅文献,不同类型的失配滤波器,并用其对LFM信号进行脉冲压缩,并观察效果。
- 去中国期刊网、IEEE网站等查找频率估计方法,利用至少两种频率估计方法编程求给定信号的频率 ,要求:
实验过程
验证DFT的对称性质
DFT 的对称性质包括:
序列 | DFT |
---|---|
$x[n]=x_\text{re}[n]+j x_\text{im}[n]$ | $X[k]=X_\text{re}[k]+j X_\text{im}[k]$ |
$x^*[n]$ | $X^*[\langle-k\rangle_N]$ |
$x^*[\langle-n\rangle_N]$ | $X^*[k]$ |
$x_\text{re}[n]$ | $X_\text{cs} =\frac{1}{2} { X[k]+X^*[\langle -k \rangle_N] }$ |
$j x_\text{im}[n]$ | $X_\text{ca}[k]=\frac{1}{2} { X[k]-X^*[\langle -k \rangle_N] }$ |
$x_\text{cs}[n]$ | $X_\text{re}[k]$ |
$x_\text{ca}[n]$ | $jX_\text{im}[k]$ |
由共轭对称性可以推导出后面四条式子,所以我们只验证后面两条式子。我们使用的 $x[n]$ 为:
\[\begin{bmatrix} 4j & 1+3j & 2+2j & 3+j & 4 \end{bmatrix}\]为了方便作图,我们定义 dftplot.m
如下:
function dftplot(xn)
%plot xn and Xk
N = length(xn);
n = 0:(N-1);
Xk = fft(xn);
subplot(321);
stem(n, real(xn));
title('real of x_n');
xlabel('n');
subplot(322);
stem(n,imag(xn));
title('image of x_n');
xlabel('n');
subplot(323);
stem(n, abs(Xk));
title('abs of X_k');
xlabel('k');
subplot(324);
stem(n, angle(Xk));
title('angle of X_k');
xlabel('k');
subplot(325);
stem(n, real(Xk));
title('real of X_k');
xlabel('k');
subplot(326);
stem(n, imag(Xk));
title('image of X_k');
xlabel('k');
end
然后我们求 $x_\text{cs}$ 和 $x_\text{ca}$,我们定义 csca.m
如下:
function [x_cs,x_ca] = csca(x)
%computes x_cs and x_ca of x
N = length(x);
n = [0:N-1];
xc = conj(x(mod(-n, N)+1));
x_cs = 0.5*(x + xc);
x_ca = 0.5*(x - xc);
end
利用这个求出 $x_\text{cs}$ 和 $x_\text{ca}$,并作出图像:
x = [4j 1+3j 2+2j 3+j 4];
[x_cs,x_ca] = csca(x);
dftplot(x);
figure;
dftplot(x_cs);
figure;
dftplot(x_ca);
$x[n]$ | $x_\text{cs}[n]$ | $x_\text{ca}[n]$ |
---|---|---|
观察最后一行可以发现,$x_\text{cs}[n]$ 的 DFT 为 $X_\text{re}[k]$,$x_\text{ca}[n]$ 的 DFT 为 $X_\text{im}[k]$。
仿照上面的过程,可以证明 $x_\text{re}[n]$ 的 DFT 为 $X_\text{cs}[k]$,$x_\text{im}[n]$ 的 DFT 为 $X_\text{ca}[k]$。先定义 idftplot.m
,方便后面作图:
subplot(321);
stem(n, real(xn));
title('real of x_n');
xlabel('n');
subplot(322);
stem(n,imag(xn));
title('image of x_n');
xlabel('n');
subplot(323);
stem(n, abs(Xk));
title('abs of X_k');
xlabel('k');
subplot(324);
stem(n, angle(Xk));
title('angle of X_k');
xlabel('k');
subplot(325);
stem(n, real(Xk));
title('real of X_k');
xlabel('k');
subplot(326);
stem(n, imag(Xk));
title('image of X_k');
xlabel('k');
然后求出 $X_\text{cs}$ 和 $X_\text{ca}$,并作出图像:
X = fft(x);
[X_cs, X_ca] = csca(X);
idftplot(X)
figure;
idftplot(X_cs);
figure;
idftplot(X_ca);
$X[k]$ | $X_\text{cs}[k]$ | $X_\text{ca}[k]$ |
---|---|---|
注意第一行,$x_\text{re}[n]$ 的 DFT 为 $X_\text{cs}[k]$,$x_\text{im}[n]$ 的 DFT 为 $X_\text{ca}[k]$。
综上,可以验证 DFT 的对称性质成立。
利用n点DFT求2n点实序列DFT
求解过程如下:
function Xk = nfft(xn)
N = length(xn);
n = 0:(N-1);
if mod(N,2) ~= 0
error('the length of xn must be 2n')
end
gn = xn(2:2:N); %偶数
hn = xn(1:2:N); %奇数
yn = gn+1i*hn; %偶数作为实部,奇数作为虚部
Yk = fft(yn); %计算 fft
[Gk,Hk] = csca(Yk); %分解为共轭对称与共轭反对称
Hk = Hk./1i;
Xk = Hk(mod(n,N/2)+1)+exp(-1i*2*pi*n/N).*Gk(mod(n,N/2)+1); %求出最后的 Xk
end
将上述代码保存为 nfft.m
,然后与 matlab 自带的 fft 进行对比:
x = [1 2 3 4 5 6];
X_fft = fft(x)
X_nfft = nfft(x)
输出的结果为:
X_fft =
1 至 3 列
21.0000 + 0.0000i -3.0000 + 5.1962i -3.0000 + 1.7321i
4 至 6 列
-3.0000 + 0.0000i -3.0000 - 1.7321i -3.0000 - 5.1962i
X_nfft =
1 至 3 列
21.0000 + 0.0000i -3.0000 + 5.1962i -3.0000 + 1.7321i
4 至 6 列
-3.0000 - 0.0000i -3.0000 - 1.7321i -3.0000 - 5.1962i
两个结果一致,说明可用 单个 N 点 DFT 计算一个实序列的 2N 点 DFT.
信号频率估计
用采样率8000Hz,对频率为350Hz的单频正弦信号进行采样,取70个点作为原始信号:
fs = 8000;
f=350;
n = 0:69;
x = sin(2*pi*f*n./fs);
加入噪声。根据搜索到的资料,可以使用 matlab 内置的函数 awgn.m
来加入高斯白噪声:
x_noise = awgn(x,20,'measured');
为了方便后续进行测试,我们定义一个测试函数(FT.m
,Frequent Test),输入频率估计方法,输出相应的均方误差:
function mse_result = FT(freq, x, snr, func,varargin)
%use func to estimate the frequent of a noise-corrupted signal, and compare it to the original freq, return mean-square-error of the estimation
% freq = original frequent
% x = original signal
% snr = sign noise ratio
% func = estimation function
% test_times = the times of test
if ~exist('test_times', 'var') || isempty(test_times)
% test_times 参数为空,或者不以变量的形式存在;
test_times=100;
end
if ~exist('fs', 'var') || isempty(fs)
% test_times 参数为空,或者不以变量的形式存在;
test_times=100;
end
errors = zeros(1,test_times);
for k = [1:test_times]
x_noise = awgn(x, snr, 'measured');
f_estimate = func(x_noise, fs);
errors(k) = [errors f_estimate-freq];
end
mse_result = sum(errors.^2) / test_times;
end
下面来讨论不同的频率估计方法。首先,最简单的,就是直接求出 DFT 最大值处对应的频率(准确来说,是 DFT 前 $N/2$ 个点中最大值),这称为谱峰法。
function f_estm = f_estm_fft(xn, fs)
Xk = abs(fft(xn));
N = length(xn);
[M,I]=max(Xk(1:ceil(N/2)));
f_estm = (N/2-abs(N/2-(I-1)))*fs/N; %matlab 编号从1开始,故要减去1
end
测试代码如下(后面的测试代码都是类似的,只需要修改 test 对应的函数即可)
snr_array = [-5 0 5 10 20 inf];
errors = [];
test = @f_estm_fft;
for snr=snr_array
errors = [errors log(FT(f, x, snr, test, 'fs',8000, 'test_times',1000))];
end
stem([-5 0 5 10 20 40], errors)
xlabel('snr/dB');
ylabel('log(mse)');
set(gca,'xtick',[-5 0 5 10 20 40]);
set(gca,'xticklabel',{'-5', '0', '5','10','20','inf'});
从下面的均方误差图可以看出,在信噪比较大的情况下,均方误差大约在 4 左右,而当信噪比较大时,均方误差就会升高:
上面的估计很简单,但存在一定误差。假设 DFT 幅度最大值的下标为 $k_p$,那么实际频率可以表示为:
\[f = \frac{k_p+\delta}{N} f_s\]我们的目标就是求出 $\delta$。常用的方法有:music算法,rife算法,周期图法,pisarenko,levenson法。下面利用 rife 算法,其表达式如下:
\[\delta = \begin{cases} \dfrac{| X(k_p+1)|}{|X(k_p+1)|+|X(k_p)|},|X(k_p+1)|>|X(k_p-1)|\\ \dfrac{-| X(k_p-1)|}{|X(k_p-1)|+|X(k_p)|},|X(k_p+1)|\leq |X(k_p-1)| \end{cases}\]function f_estm = f_estm_rife(xn, fs)
Xk = abs(fft(xn));
N = length(xn);
[M,k0]=max(Xk(1:ceil(N/2)));
p0=abs(Xk(k0));
p1=abs(Xk(k0-1));
p2=abs(Xk(k0+1));
if p1>=p2
pp=p1;flag=-1;
else
pp=p2;flag=1;
end
k=k0-1+flag*pp/(p0+pp);
f_estm=k*fs/N;
end
使用相同的测试代码得到的图像如下,可以看出在信噪比越大,其估计的准确度越好。
另外还使用课堂中提到的导数法,其公式如下:
\[\hat{\omega_c} = \frac{\sum_{k=0}^{N-1}|X_k|^2 k\omega_0}{\sum_{k=0}^{N-1}|X_k|^2}\approx \frac{\sum_{k=k_0-\Delta}^{k_0+\Delta}|X_k|^2 k\omega_0}{\sum_{k=k_0-\Delta}^{k_0+\Delta}|X_k|^2}\]function f_estm = f_estm_dif(xn, fs)
Xk = abs(fft(xn));
N = length(xn);
[M,k0]=max(Xk(1:ceil(N/2)));
Delta=4;
num=0;
den=0;
for k = [k0-Delta:k0+Delta]
if k<1 || k>N
continue
end
num=num+abs(Xk(k))^2.*(k-1);
den=den+abs(Xk(k))^2;
end
f_estm=num*fs/(den*N);
end
使用相同的测试代码得到的图像如下,可以看出导数法在噪声较大与 rife 法类似,而在噪声较小时不如 rife 法,总体上仅比谱峰法好一点点。
总结:snr≤5时都不好,-5<snr<10 时谱峰法好,snr>10时rife法好啊👍
LFM脉冲压缩
说真的,我压根就不知道这是什么玩意,也没时间慢慢推导了,就直接参考了这篇文章:线性调频(LFM信号)脉冲压缩雷达matlab仿真,以及 基于MATLAB的线性调频信号的仿真
先是生成线性调频信号并加入噪声:
%%线性调频信号
T=10e-6; %p脉冲持续时间10us
B=30e6; %线性调频信号的频带宽度30MHz
K=B/T; %调频斜率
Fs=2*B;Ts=1/Fs; %采样频率和采样间隔
N=T/Ts;
t=linspace(-T/2,T/2,N);
St=exp(j*pi*K*t.^2); %线性调频信号
St=awgn(St,10);
subplot(211)
plot(t*1e6,real(St));
xlabel('时间/us');
title('线性调频信号的实部');
grid on;axis tight;
subplot(212)
freq=linspace(-Fs/2,Fs/2,N);
plot(freq*1e-6,fftshift(abs(fft(St))));
xlabel('频率/MHz');
title('线性调频信号的幅频特性');
grid on;axis tight;
%% 线性调频与脉冲压缩
clear,clc,close all
set(0,'defaultfigurecolor','w')
%% Chirp信号参数设置
Tr = 10e-6;%时宽
Br = 30e6;%带宽
Fs = 2*Br;%采样率
%% Chirp信号参数导出
Kr = Br/Tr;%调频率
N = round( Tr / (1/Fs) );%采样点数
t = linspace( -Tr/2+0.1*Tr , Tr/2+0.1*Tr , N);%在[-Tp/2,Tp/2]选取采样点
%% Chirp信号生成
st = ( abs(t) < Tr/2 ) .* exp( 1j * pi * Kr * t.^2 );
st=awgn(st,10);
f_chirp= Kr * t; %信号频率
phase_chirp = pi * Kr * t.^2;%信号相位
%% 频谱
freq = linspace(-Fs/2,Fs/2,N);%频域采样
Sf = fftshift( fft(st) );
%% 时域匹配滤波
ht = conj( fliplr(st) ); %时域匹配滤波为发射信号时间反褶再取共轭
s1 = conv(st,ht); %线性调频信号经过匹配滤波器后的输出(时域卷积)
N1 = N+N-1 ;%线性卷积后信号长度变为 N1+N2-1
t1 = linspace( -Tr/2 , Tr/2 , N1);
% 时域匹配滤波
figure,plot( t1*1e6 , abs(s1) ),xlabel('t /us'),ylabel('幅度谱'),title('时间反褶取共轭,时域卷积');
%%脉冲压缩
clear;
close all;
T=10e-6; %时宽
B=30e6; %带宽
Rmin=8500;
Rmax=11500;
R=[9000,10000,10200];
RCS=[1 1 1];
C=3e8;
K=B/T;
Rwid=Rmax-Rmin;
Twid=2*Rwid/C;
Fs=5*B;Ts=1/Fs;
Nwid=ceil(Twid/Ts);
t=linspace(2*Rmin/C,2*Rmax/C,Nwid);
M=length(R);
td=ones(M,1)*t-2*R'/C*ones(1,Nwid);
SNR=[1,0.001];
for i=1:1:2
Srt1=RCS*(exp(1i*pi*K*td.^2).*(abs(td)<T/2));
n=sqrt(0.5*SNR(i))*(randn(size(Srt1))+1i*randn(size(Srt1))); %加噪声
Srt=Srt1+n;
Nchirp=ceil(T/Ts);
Nfft=2^nextpow2(Nwid+Nwid-1);
Srw=fft(Srt,Nfft);
rw1=fft(Srt1,Nfft);
t0=linspace(-T/2,T/2,Nchirp);
St=exp(1i*pi*K*t0.^2);
Sw=fft(St,Nfft);
Sot=fftshift(ifft(Srw.*conj(Sw)));
N0=Nfft/2-Nchirp/2;
Z=abs(Sot(N0:N0+Nwid-1));
figure
subplot(211)
plot(t*1e6,real(Srt));
axis tight;
xlabel('us');
ylabel('幅度')
title(['加噪线性调频信号压缩前,SNR =',num2str(-1*10*log10(SNR(i)))]);
subplot(212)
plot(t*C/2,Z)
xlabel('Range in meters');ylabel('幅度 ')
title(['加噪线性调频信号压缩后,SNR =',num2str(-1*10*log10(SNR(i)))]);
end
最后这个实验到底是什么鬼啊😥