KerberosSDR代码笔记(5) 信号处理（采样时间延迟计算、相位差计算的2种方法、MUSIC算法）

https://github.com/rtlsdrblog/kerberossdr/blob/master/_signalProcessing/hydra_signal_processor.py

这部分是最重要的部分，它里面包含几个重要的算法，采样时间同步，相位同步，以及调用了空间谱估计算法。

# KerberosSDR Signal Processor
#
# Copyright (C) 2018-2019  Carl Laufer, Tamás Pető
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
#
#
# -*
# - coding: utf-8 -*-import sys
import os
import time
# Math support
import numpy as np# Signal processing support
from scipy import fft,ifft
from scipy import signal
from scipy.signal import correlate #互相关函数# Plot support
#import matplotlib.pyplot as plt# GUI support
from PyQt4 import QtGui, QtCore# Import the pyArgus module
#root_path = os.getcwd()
#pyargus_path = os.path.join(os.path.join(root_path, "pyArgus"), "pyArgus")
#sys.path.insert(0, pyargus_path)
#import directionEstimation_v1p15 as defrom pyargus import directionEstimation as de #测向算法# Import APRiL module
#april_path = os.path.join(os.path.join(root_path, "_APRIL"), "APRIL")
#sys.path.insert(0, april_path)
#import channelPreparation as cp
#import clutterCancellation as cc
#import detector as detfrom pyapril import channelPreparation as cp #这些都是被动雷达的
from pyapril import clutterCancellation as cc
from pyapril import detector as det
from pyapril.hitProcessor import CA_CFARclass SignalProcessor(QtCore.QThread):#初始化了几个信号signal_spectrum_ready = QtCore.pyqtSignal() #频谱计算完成的信号signal_sync_ready = QtCore.pyqtSignal() #采样延迟计算完成的信号signal_DOA_ready = QtCore.pyqtSignal() #测向计算完成的信号signal_overdrive = QtCore.pyqtSignal(int) #信号是否饱和的信号，但是没用到signal_period    = QtCore.pyqtSignal(float) #计算周期signal_PR_ready = QtCore.pyqtSignal()def __init__(self, parent=None, module_receiver=None):"""Description:------------Parameters:-----------Return values:--------------"""super(SignalProcessor, self).__init__(parent)self.module_receiver = module_receiver #传入接收机对象给新处理模块self.en_spectrum = True  #是否显示频谱self.en_sync = True  #这里对应的是界面上的显示同步曲线（enable sync display）self.en_sample_offset_sync = False#这里对应的是界面上的执行同步按钮（sample sync)self.en_record = False #是否保存数据self.en_calib_iq = False #是否做相位差校准self.en_calib_DOA_90 = Falseself.en_DOA_estimation = False #是否要做测向self.en_PR_processing = Falseself.en_PR_autodet = False# DOA processing optionsself.en_DOA_Bartlett = Falseself.en_DOA_Capon = Falseself.en_DOA_MEM = Falseself.en_DOA_MUSIC = False #这是选择使用哪个测向算法，一般用MUSICself.en_DOA_FB_avg = False #这是是否打开FB averageself.DOA_inter_elem_space = 0.5 #这是lamdaself.DOA_ant_alignment = "ULA" #这是阵型# Passive Radar processing parameters #这都是被动雷达用的参数self.ref_ch_id = 0self.surv_ch_id = 1self.en_td_filtering = Falseself.td_filter_dimension = 1        self.max_Doppler = 500  # [Hz]self.windowing_mode = 0self.max_range = 128  # [range cell]self.cfar_win_params = [10,10,4,4] # [Est. win length, Est. win width, Guard win length, Guard win width]self.cfar_threshold = 13self.RD_matrix = np.ones((10,10))self.hit_matrix = np.ones((10,10))self.RD_matrix_last = np.ones((10,10))self.RD_matrix_last_2 = np.ones((10,10))self.RD_matrix_last_3 = np.ones((10,10))self.center_freq = 0  # TODO: Initialize this [Hz] #中心频率self.fs = 1.024 * 10**6  # Decimated sampling frequncy - Update from GUI#这是经过降采样的采样率#self.sample_size = 2**15self.channel_number = 4 #接收机通道数量# Processing parameters        self.test = Noneself.spectrum_sample_size = 2**14 #2**14 #用于频谱运算的采样点数量self.DOA_sample_size = 2**15 # Connect to GUI value?? #用于测向运算的采样点数量self.xcorr_sample_size = 2**18 #2**18 #用于互相关运算的采样点数量self.spectrum = np.ones((self.channel_number+1,self.spectrum_sample_size), dtype=np.float32)   #这里存储的是频谱的数据，每个通道一行，所以行数是channel_number，#但是第0行要存储横轴的频率值信息，所以还要+1#列的数量就是频谱显示的采样点数量self.xcorr = np.ones((self.channel_number-1,self.xcorr_sample_size*2), dtype=np.complex64) #这里存储的是互相关数据，都是另3个通道与第1个通道互相关运算的结果#这个结果只有3个，所以只有channel_number-1行，互相关是卷积，长度会补得比函数本身长#所以是原函数的采样点数量*2 xcorr_sample_size*2self.phasor_win = 2**10 # Phasor plot window#这个本来用来绘制相位差的但是实际没用self.phasors = np.ones((self.channel_number-1, self.phasor_win), dtype=np.complex64)self.run_processing = False #这是用来控制主循环的标志位# Result vectorsself.delay_log= np.array([[0],[0],[0]]) #记录采样时间延迟的数组，另3个通道相对第1个通道的，每个结果放在一行里，所以有3行#每一列是不同时间点求的延迟结果self.phase_log= np.array([[0],[0],[0]])#记录相位差的数组，类似采样时间延迟，也是另3个通道相对第一个通道的，有3行#每一列是不同时间点求的相位差的结果self.DOA_Bartlett_res = np.ones(181)self.DOA_Capon_res = np.ones(181)self.DOA_MEM_res = np.ones(181)self.DOA_MUSIC_res = np.ones(181)#这里存储的是MUSIC算法计算出的功率谱密度函数self.DOA_theta = np.arange(0,181,1)#这里存储的是DOA测出的角度# Auto resync params #自动定期同步用的参数，实际没有实现self.lastTime = 0self.runningSync = 0self.timed_sync = Falseself.noise_checked = Falseself.resync_time = -1def run(self):# xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx#    # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxself.run_processing = Truewhile self.run_processing:start_time = time.time() #记录进入循环的时间点# Download samples#if(self.en_sync or self.en_spectrum): time.sleep(0.25) # You can play with this value, but it may affect stability   self.module_receiver.download_iq_samples()#告诉receiver程序，可以下载iq数据了，它会通知c程序开门#这样4个接收机通道的数据都会存进iq_samples这个二维数组里self.DOA_sample_size = self.module_receiver.iq_samples[0,:].size#用iq_samples中得到的实际长度（第1个接收机通道的数据长度）来决定DOA计算的长度self.xcorr_sample_size = self.module_receiver.iq_samples[0,:].size#互相关计算的长度也是iq_samples的第1个通道的实际长度self.xcorr = np.ones((self.channel_number-1,self.xcorr_sample_size*2), dtype=np.complex64) #这里是存储互相关结果的地方，前面说了，互相关结果是后3个通道与第1个通道做的#所以能得到3组结果，就是3行，所以行数是channel_number-1#但是互相关结果的长度由于是卷积运算，所以比原来的函数长，要xcorr_sample_size*2#原函数里的采样点都是复数，互相关结果也是复数，所以dtype是np.complex64# Check overdrive #检测信号饱和，实际没有用if self.module_receiver.overdrive_detect_flag:self.signal_overdrive.emit(1)else:self.signal_overdrive.emit(0)# Display spectrum #显示频谱if self.en_spectrum:  #先判断是否勾选了显示频谱self.spectrum[0, :] = np.fft.fftshift(np.fft.fftfreq(self.spectrum_sample_size, 1/self.fs))/10**6#self.spectrum的第0行存储的应该是频谱的频率范围#np.fft.fftfreq可以生成这个范围，它的参数是fft点数（输入和输出的数据长度）#和采样周期（也就是采样率的倒数）#np.fft.fftshift的作用是将0频点移到中间#最后还要除以10**6，只保存MHz前的数字m = self.channel_number#self.spectrum[1:m+1,:] = 10*np.log10(np.fft.fftshift(np.abs(np.fft.fft(self.module_receiver.iq_samples[0:m, 0:self.spectrum_sample_size]))))for m in range(self.channel_number):self.spectrum[m+1,:] = 10*np.log10(np.fft.fftshift(np.abs(np.fft.fft(self.module_receiver.iq_samples[m, 0:self.spectrum_sample_size]))))#用一个循环来生成4个接收机通道的fft结果，并存入self.spectrum的后几行#首先读出iq_samples里第m行的数据，代表第m个接收机的采样点#纵轴从0开始表示第一个数据开始算#直到达到要用来算频谱的数据长度spectrum_sample_size#有了数据后就用np.fft.fft来计算fft结果#再求np.abs这样只保留实数值，就是我们熟悉的频谱图了#然后再用np.fft.fftshift把零频点移到频谱中间#再用10*np.log10把单位转为dBself.signal_spectrum_ready.emit() #发出信号表示频谱算好了# Synchronization #同步if self.en_sync or self.timed_sync: #print("Sync graph enabled") #这里勾选了同步显示self.sample_delay() #一旦勾选同步显示，就会计算几个通道的延迟大小self.signal_sync_ready.emit() #上面函数返回表示计算完成了，发出信号# Sample offset compensation request #延迟补偿if self.en_sample_offset_sync: #这里是点击了同步按钮self.module_receiver.set_sample_offsets(self.delay_log[:,-1])#把之前计算到的各通道的延迟大小的最新值应用到接收机python程序里#[:,-1]代表所有行的最后一列数据self.en_sample_offset_sync = False #纠正完了采样延迟就可以把这个按钮取消了# IQ calibration request #相位差校准if self.en_calib_iq: #点击相位差校准按钮# IQ correctionfor m in range(self.channel_number):self.module_receiver.iq_corrections[m] *= np.size(self.module_receiver.iq_samples[0, :])/(np.dot(self.module_receiver.iq_samples[m, :],self.module_receiver.iq_samples[0, :].conj()))#这里在计算每个接收机通道需要的校准值，这个校准值的分母是np.dot两个复数点乘#2个复数点乘得到的复数的相位是这2个复数的相位之和#第一个复数是m通道的采样点#第二个复数是第1个接收机通道的共轭（相当于相位加了负号）#这样点乘得到的复数的相位就是m通道相对第1个接收机通道的相位差#把含有相位差信息的复数放到分母上，这样其它通道乘以它就会把自己相位减去相位差#分子是第一个通道的采样点大小，应该是用来归一化用的？#iq_correction原本的值是1，所以第一次运算*=就跟=直接赋值是一样的#后期再做这一步就是把最新的相位差继续乘上去#相当于继续纠正相位差（但一般一次就够了，不需要后续再做相位差校准）c = np.sqrt(np.sum(np.abs(self.module_receiver.iq_corrections)**2))#这里先对前面算出来的iq_correction求模，再求平方，又把数组全部加起来后开根号#得到的c应该也是归一化用的？self.module_receiver.iq_corrections = np.divide(self.module_receiver.iq_corrections, c)#用np.divide把之前求的iq_corrections除以c就得到了要用的相位差校准值#print("Corrections: ",self.module_receiver.iq_corrections)self.en_calib_iq = False #得到相位差校准值后就可以把相位差校准按钮取消了if self.en_calib_DOA_90: #这个暂时还在实验阶段#TODO: Experimental only for UCA, implement this properly!# This calibration is currently done for 0 deg not 90 x = self.DOA_inter_elem_space * np.cos(2*np.pi/4 * np.arange(4))y = self.DOA_inter_elem_space * np.sin(-2*np.pi/4 * np.arange(4)) # For this specific array onlyref_vector = de.gen_scanning_vectors(4, x, y, np.zeros(1))[:, 0]                #ref_vector = np.exp(1j*2*np.pi*0.5*np.cos(np.radians(0-np.arange(self.channel_number)*(360)/self.channel_number))) # UCA                N= np.size(self.module_receiver.iq_samples[0, :])for m in range(self.channel_number):self.module_receiver.iq_corrections[m] *= ref_vector[m]*N/(np.dot(self.module_receiver.iq_samples[m, :],self.module_receiver.iq_samples[0, :].conj()))                #print("Corrections: ",self.module_receiver.iq_corrections)self.en_calib_DOA_90 = False# Direction of Arrival estimation #DOA计算if self.en_DOA_estimation: #如果勾选了DOA计算# Get FFT for squelchself.spectrum[1,:] = 10*np.log10(np.fft.fftshift(np.abs(np.fft.fft(self.module_receiver.iq_samples[0, 0:self.spectrum_sample_size]))))#获取第1个接收机通道的fft，用来滤除小信号self.estimate_DOA() #计算DOAself.signal_DOA_ready.emit() #计算完DOA要发出信号# Passive Radar processing #被动雷达处理if self.en_PR_processing:
#                self.module_receiver.channel_number = 2self.PR_processing()self.signal_PR_ready.emit()
#            else:
#                self.module_receiver.channel_number = 4# Record IQ samples #保存iq采样点if self.en_record: #如果打开了会把iq_samples里的采样点存下来np.save('hydra_samples.npy', self.module_receiver.iq_samples)# Code to maintain sync #定期同步，但是没有使用'''if self.timed_sync and not self.en_sync:if not self.noise_checked:self.module_receiver.switch_noise_source(0)self.timed_sync = Falseself.en_sample_offset_sync=Trueself.runningSync = 0resync_on = Trueif(self.resync_time < 10):resync_on = Falseif(((start_time - self.lastTime) > self.resync_time) and not self.en_sync and resync_on):self.lastTime = start_timeself.module_receiver.switch_noise_source(1)time.sleep(0.1)self.runningSync = 1self.timed_sync = True'''stop_time = time.time() #循环结束前计时self.signal_period.emit(stop_time - start_time) #这里会发出信号，并且包含测向运行时所有的获取数据+运算需要的时间多长def sample_delay(self): #计算采样时间延迟#print("Entered sample delay func")N = self.xcorr_sample_size #互相关大小iq_samples = self.module_receiver.iq_samples[:, 0:N] #根据要计算的互相关长度，从iq_samples里取出4个接收机通道的数据#互相关结果是xcorr_sample_size的两倍，xcorr_sample_size只是互相关的输入原函数长度#xcorr_sample_size=2**18=256*1024，iq_samples里的数据正好够用delays = np.array([[0],[0],[0]]) #延迟phases = np.array([[0],[0],[0]]) #相位差#上面的2个数组看着都是二维，实际都是一维数组，是3行1列的数组#delays只存储3个值，是目前的另3个通道相对第1通道的延迟#都是当前值，不包含历史信息，phases也是如此# Channel matchingnp_zeros = np.zeros(N, dtype=np.complex64)#用来填补互相关原函数用的#互相关运算是2个原函数求卷积，卷积结果的长度是原函数长度的2倍，原函数不够长#需要用0把它们填满，使其长度跟结果一样长，所以np_zeros长度也是Nx_padd = np.concatenate([iq_samples[0, :], np_zeros])#这个concatenate就是用来把原函数和np_zeros接起来的#根据卷积运算的要求，两个原函数一个从左往右，另一个从右往左，各个点相乘再求和#所以x_padd是把第1个通道的值取出来再把np_zeros填在右侧，并且接起来x_fft = np.fft.fft(x_padd)#为了计算方便，在原始数据数据处理完后选择在频域计算#频域只要求一个原函数乘以另一个原函数的共轭就得到互相关的fft，再用ifft变回时域就行#所以这里用x_padd做了fft得到了x_fftfor m in np.arange(1, self.channel_number):y_padd = np.concatenate([np_zeros, iq_samples[m, :]])y_fft = np.fft.fft(y_padd)#类似的，这里对第m通道的采样数据求了y_padd，然后再求了y_fftself.xcorr[m-1] = np.fft.ifft(x_fft.conj() * y_fft)#这样就得到了第m通道相对于第1个接收机通道(m=0)的互相关函数delay = np.argmax(np.abs(self.xcorr[m-1])) - N#互相关函数的模的最大值对应横坐标就是延迟大小，为什么要减去N，暂时不知道？#这个可能是类似fftshift类似，延迟可能是正可能是负（提前）#但是目前的互相关函数的中间点应该不在0点而是在N点，所以要减去N#相当于把函数往左搬回0点，得到的delay才是真正的延迟#phase = np.rad2deg(np.angle(self.xcorr[m-1, delay + N]))phase = np.rad2deg(np.angle(self.xcorr[m-1, N]))#互相关函数第N个位置，也就是中间点或者理论上对应0位置的点的复数角就是相位差？#这里还不是很清楚#offset = 50000                     #self.phasors[m-1, :] = (iq_samples[0, offset: self.phasor_win+offset] * iq_samples[m, offset+delay: self.phasor_win+offset+delay].conj())#self.phasors[m-1, :] = (iq_samples[0, 0: self.phasor_win] * iq_samples[m, 0: self.phasor_win].conj())#本来记录用的相位差也是把一个采样点的值乘以另一个采样点的共轭求出来的#与相位差校准时原理一样，但是现在不用了，直接用互相关结果的复数来得到？"""self.IQSamples[1, :] = np.roll(self.IQSamples[1, :], delay * -1)if delay > 0:self.IQSamples[1, -delay::] = np.zeros(delay, dtype=np.complex64)if delay < 0:self.IQSamples[1, 0: np.abs(delay)] = np.zeros(np.abs(delay), dtype=np.complex64)"""#msg = "[ INFO ] delay: " + str(delay)#print(msg)delays[m-1,0] = delay phases[m-1,0] = phase#计算完成后就要把计算结果存入3行1列的数组了，列号都是0，行号是当前计算通道-1#比如第2个接收机通道(m=1)相对于第1个通道的延迟存入delays[0,0]是第一行第一列#第3个接收机通道(m=2)存入delays[1,0]第二行第一列#第4个接收机通道(m=3)存入delays[2,0]第三行第一列#相位差phases也是类似self.delay_log = np.concatenate((self.delay_log, delays),axis=1)self.phase_log = np.concatenate((self.phase_log, phases),axis=1)#得到当前时间点的数据后，就要把它们存入历史记录了#直接用np.concatenate把2个3行1列的数组接到delay_log和phase_log的最后就行了def delete_sync_history(self): #删除同步历史信息self.delay_log= np.array([[0],[0],[0]]) #清空采样时间延迟数组self.phase_log= np.array([[0],[0],[0]]) #清空相位差数组#上面2个都是二维数组，有3行，每行好多列，分别代表另3个通道相对第1个通道的#采样时间延迟和相位差数据，有好多列是因为把不同时间算出的数据都存下来了#越靠后的数据越新def estimate_DOA(self): #计算DOA，主要都在调用pyargus的de包#print("[ INFO ] Python DSP: Estimating DOA")iq_samples = self.module_receiver.iq_samples[:, 0:self.DOA_sample_size]#首先要截取iq_samples，4个通道的数据都要，所以每行都要了，行用了冒号表示#但是DOA算法只需要计算DOA_sample_size长度的数据，所以列是0:self.DOA_sample_size# Calculating spatial correlation matrix #计算协方差矩阵R = de.corr_matrix_estimate(iq_samples.T, imp="fast")if self.en_DOA_FB_avg: #如果是线阵并且打开了FB average要更新协方差矩阵R=de.forward_backward_avg(R)M = np.size(iq_samples, 0) #iq_samples的行数，就是接收机通道数，4if self.DOA_ant_alignment == "UCA": #均匀圆阵self.DOA_theta =  np.linspace(0,360,361)#scanning_vectors = de.gen_uca_scanning_vectors(M, self.DOA_inter_elem_space, self.DOA_theta)x = self.DOA_inter_elem_space * np.cos(2*np.pi/M * np.arange(M))y = self.DOA_inter_elem_space * np.sin(-2*np.pi/M * np.arange(M)) # For this specific array onlyscanning_vectors = de.gen_scanning_vectors(M, x, y, self.DOA_theta)#生成扫描向量# DOA estimationif self.en_DOA_Bartlett:self.DOA_Bartlett_res = de.DOA_Bartlett(R, scanning_vectors)if self.en_DOA_Capon:self.DOA_Capon_res = de.DOA_Capon(R, scanning_vectors)if self.en_DOA_MEM:self.DOA_MEM_res = de.DOA_MEM(R, scanning_vectors,  column_select = 0)if self.en_DOA_MUSIC:self.DOA_MUSIC_res = de.DOA_MUSIC(R, scanning_vectors, signal_dimension = 1)#使用MUSIC算法，计算协方差矩阵和扫描向量，得到DOA结果elif self.DOA_ant_alignment == "ULA": #均匀线阵self.DOA_theta =  np.linspace(-90,90,181)x = np.zeros(M)y = np.arange(M) * self.DOA_inter_elem_space            scanning_vectors = de.gen_scanning_vectors(M, x, y, self.DOA_theta)# DOA estimationif self.en_DOA_Bartlett:self.DOA_Bartlett_res = de.DOA_Bartlett(R, scanning_vectors)if self.en_DOA_Capon:self.DOA_Capon_res = de.DOA_Capon(R, scanning_vectors)if self.en_DOA_MEM:self.DOA_MEM_res = de.DOA_MEM(R, scanning_vectors,  column_select = 0)if self.en_DOA_MUSIC:self.DOA_MUSIC_res = de.DOA_MUSIC(R, scanning_vectors, signal_dimension = 1)#print(self.DOA_MUSIC_res)def PR_processing(self): #被动雷达的算法#print("[ INFO ] Python DSP: Start Passive Radar processing")ref_ch = self.module_receiver.iq_samples[self.ref_ch_id, :]surv_ch = self.module_receiver.iq_samples[self.surv_ch_id, :]if self.en_td_filtering:surv_ch, w = cc.Wiener_SMI_MRE(ref_ch, surv_ch, self.td_filter_dimension)
#            surv_ch, w = cc.fast_wiener(self.td_filter_dimension, ref_ch, surv_ch)#surv_ch, w = cc.Wiener_SMI(ref_ch, surv_ch, self.td_filter_dimension, imp="fast")#print("[ DONE ] Timde domain filtering finished")if(self.windowing_mode == 0):pass#surv_ch = det.windowing(surv_ch, "Rectangular")else:surv_ch = det.windowing(surv_ch, "Hamming")self.RD_matrix = det.cc_detector_ons(ref_ch, surv_ch, self.fs, self.max_Doppler, self.max_range, verbose=0, Qt_obj=None)if self.en_PR_autodet:self.hit_matrix = CA_CFAR(self.RD_matrix,self.cfar_win_params, self.cfar_threshold)            #print("[ DONE ] Range-Doppler processing finished")def stop(self):self.run_processing = Falsedef busy_wait(dt):current_time = time.time()while (time.time() < current_time+dt):pass

接下来看一下最常用的UCA均匀圆阵所涉及的几个函数的具体实现：

de.corr_matrix_estimate，de.gen_scanning_vectors，de.DOA_MUSIC

它们都来自于pyargus库

https://github.com/petotamas/pyArgus

这是一个跟天线阵有关的库，有3个主要的文件，antennaArrayPattern.py，这是用来绘制阵列天线方向图的，beamform.py，用来实现波束成形，directionEstimation.py，计算信号到达角的，除此之外还有一些用来测试的程序以及配套文档。

这几个模块互相之间没有依赖，所以我们可以直接看directionEstimation.py

https://github.com/petotamas/pyArgus/blob/master/pyArgus/directionEstimation.py

以及它的文档：

https://github.com/petotamas/pyArgus/blob/master/docs/nb_direction_of_arrival_estimaton.ipynb

ipynb可以用github直接查看，但是网速不好可能打不开，可以尝试这个网站nbviewer.jupyter.org

测向原理可以参考这篇文章：

https://blog.csdn.net/qq_23947237/article/details/82318222

它解释了为什么经过它提出的步骤就能测出方向。

这个步骤总结一下就是：

1.先求出阵列矩阵(scanning vector)

阵列矩阵跟阵形有关，上面那个博客里用的是线阵，但是我们实际python代码里用的是更一般化的代码，只需要填写x,y的信息，可以生成面阵的阵列矩阵

2.求协方差矩阵(correlation matrix)

可以用接收机通道上收到的大量的采样点来估计出这个协方差矩阵

3.求出噪声子空间

对协方差矩阵做特征值分解，得到这个噪声子空间

4.计算空间谱函数

利用阵列矩阵和噪声子空间计算空间谱函数，就能得到music的波形，最后做谱峰搜索就能得到入射角（程序其它部分做的）。

下面是我从hydra_signal_processing.py里摘录的使用UCA和MUSIC算法时会调用到的代码，基本就是按照上面的步骤做的，只不过把第一步和第二步换了一下顺序，这个没有影响。

iq_samples = self.module_receiver.iq_samples[:, 0:self.DOA_sample_size]
#首先要截取iq_samples，4个通道的数据都要，所以每行都要了，行用了冒号表示
#但是DOA算法只需要计算DOA_sample_size长度的数据，所以列是0:self.DOA_sample_sizeR = de.corr_matrix_estimate(iq_samples.T, imp="fast")
#使用大量的采样点来估计协方差矩阵M = np.size(iq_samples, 0)
#iq_samples的行数，就是接收机通道数，4
self.DOA_theta =  np.linspace(0,360,361)
#生成0~360，一共361个数字的等间隔数列，self.DOA_theta=[0 1 2 3 ... 360]
x = self.DOA_inter_elem_space * np.cos(2*np.pi/M * np.arange(M))
y = self.DOA_inter_elem_space * np.sin(-2*np.pi/M * np.arange(M))
#DOA_inter_elem_space里是lamda 天线间距/波长 对于圆阵还要除以sqrt(2)，晚点确认？
#np.arange(4)会生成一个数组[0 1 2 3]
#如果暂时忽略掉self.DOA_inter_elem_space
#只考虑np.cos(2*np.pi/M * np.arange(M))和np.sin(-2*np.pi/M * np.arange(M))
#它们会生成单位圆与横轴右侧交点开始逆时针旋转，均匀分布的4个点，x和y分别是这些点的横坐标和纵坐标
#[x,y] = [1,0] [0,-1] [-1,0] [0,1] 实际x和y还要乘以lamda
scanning_vectors = de.gen_scanning_vectors(M, x, y, self.DOA_theta)
#根据阵列的空间布局（包括阵形和lamda)来计算阵列向量self.DOA_MUSIC_res = de.DOA_MUSIC(R, scanning_vectors, signal_dimension = 1)
#根据协方差矩阵和阵列向量来做MUSIC计算
#signal_dimension是待测信号数量，这里规定了只有1个待测信号

可以看到hydra_signal_processing.py在做music算法时虽然步骤和理论差不多，但是只是调用了几个pyargus里的函数，具体的操作在pyargus中，我摘录了对应的几个函数。

def corr_matrix_estimate(X, imp="mem_eff"):"""Estimates the spatial correlation matrix with sample averaging    Implementation notes:--------------------Two different implementation exist for this function call. One of them use a for loop to iterate through thesignal samples while the other use a direct matrix product from numpy. The latter consumes more memory(as all the received coherent multichannel samples must be available at the same time)but much faster for large arrays. The implementation can be selected using the "imp" function parameter.Set imp="mem_eff" to use the memory efficient implementation with a for loop or set to "fast" in order to usethe faster direct matrix product implementation.Parameters:-----------:param X : Received multichannel signal matrix from the antenna array.         :param imp: Selects the implementation method. Valid values are "mem_eff" and "fast". The default value is "mem_eff".:type X: N x M complex numpy array N is the number of samples, M is the number of antenna elements.:type imp: stringReturn values:-------------:return R : Estimated spatial correlation matrix:rtype R: M x M complex numpy array:return -1 : When unidentified implementation method was specified"""      #由于调用这个函数的时候，iq_sample做了转置(iq_sample.T)#因此X变量里现在的列数是4列，对应4个接收机通道，行数是每一个接收机通道的采样点数量N = np.size(X, 0) #每个通道采样点数量M = np.size(X, 1) #接收机通道数，4个R = np.zeros((M, M), dtype=complex)    # --input check--if N < M:print("WARNING: Number of antenna elements is greather than the number of time samples")  #如果采样点比接收机通道少，那肯定不行print("WARNING: You may flipped the input matrix")# --calculation--if imp == "mem_eff":            for n in range(N):R += np.outer(X[n, :], np.conjugate(X[n, :]))elif imp == "fast":X = X.T #现在X矩阵又做了一次转置，那么X里面现在就是原始的iq_samples了#根据协方差矩阵的估计公式（参考那篇博客的文章）#接收机采样点矩阵乘以自己的共轭转置，再累加，再除以采样点数量就行#矩阵右上角的H就代表共轭转置#下面这一行做的就是采样点乘以自己的共轭转置并累加了R = np.dot(X, X.conj().T)else:print("ERROR: Unidentified implementation method")print("ERROR: No output is generated")return -1#下面这一行做的是除以采样点数量，做完了就得到协方差矩阵的估计值了R = np.divide(R, N)return Rdef gen_scanning_vectors(M, x, y, thetas):"""Description:------------This function prepares scanning vectorors for general antenna array configurations        #这里生成的是广义的阵列天线的阵列矩阵，更接近面阵或者任意平面阵Parameters:-----------:param M : Number of antenna elements on the circle:param x : x coordinates of the antenna elements on a plane:param y : y coordinates of the antenna elements on a plane:param thetas : A vector containing the incident angles e.g.: [0deg, 1deg, 2deg, ..., 180 deg]:type M: int:type x: 1D numpy array:type y: 1D numpy array:type R: float:type thetas: 1D numpy arrayReturn values:-------------:return scanning_vectors : Estimated signal dimension:rtype scanning_vectors: 2D numpy array with size: M x P, where P is the number of incident angles"""scanning_vectors = np.zeros((M, np.size(thetas)), dtype=complex)#np.zeros会生成4行，361列的数组for i in range(np.size(thetas)):        scanning_vectors[:,i] = np.exp(1j*2*np.pi* (x*np.cos(np.deg2rad(thetas[i])) + y*np.sin(np.deg2rad(thetas[i]))))    #我暂时还没找到面阵的阵列矩阵的公式，但是那篇博客里的线阵公式看着也差不多#也是e^(-j*(N)*2*pi*d*sin(theta)/lamda)，复数表示差不多，但是公式里lamda在分母#我这里x,y里包含了lamda，但是它们在分子上return scanning_vectorsdef DOA_MUSIC(R, scanning_vectors, signal_dimension, angle_resolution = 1):"""                                 MUSIC - Multiple Signal Classification methodDescription:------------    The function implements the MUSIC method for direction estimationCalculation method : 1ADORT(theta) = ---------------------------H        H S(theta) * En En  * S(theta)#这个公式和那篇博客文章里的空间谱函数一样，现在S(theta)也就是阵列矩阵已经直接传入了#我们只需要算出En就行，它是噪声子空间Parameters:-----------                :param R: spatial correlation matrix            :param scanning_vectors : Generated using the array alignment and the incident angles                                                 :param signal_dimension:  Number of signal sources    :type R: 2D numpy array with size of M x M, where M is the number of antennas in the antenna system                            :tpye scanning vectors: 2D numpy array with size: M x P, where P is the number of incident angles                                   :type signal_dimension: intReturn values:--------------:return  ADORT : Angular dependent orthogonality. Expresses the orthongonality of the current steering vector to the noise subspace:rtype : numpy array:return -1, -1: Input spatial correlation matrix is not quadratic:return -2, -2: dimension of R not equal with dimension of the antenna array :return -3, -3: Spatial correlation matrix is singular"""# --- Parameters ---  # --> Input checkif np.size(R, 0) != np.size(R, 1): #协方差矩阵应该是一个方阵print("ERROR: Correlation matrix is not quadratic")return -1, -1if np.size(R, 0) != np.size(scanning_vectors, 0): #阵列矩阵的列数要和协方差矩阵相同，都是通道数4print("ERROR: Correlation matrix dimension does not match with the antenna array dimension")return -2, -2                    ADORT = np.zeros(np.size(scanning_vectors, 1),dtype=complex)#这是计算出的谱函数，横坐标就是所有角度0~360，就是阵列矩阵的列数M = np.size(R, 0) #用协方差矩阵的行数得到通道数# --- Calculation ---# Determine eigenvectors and eigenvaluessigmai, vi = lin.eig(R) #求协方差矩阵R的特征值sigmai和特征向量vi#sigmai 是一行数组，每一元素是一个特征值#vi 也是一行，但是是数组套数组，每一个元素都是一个特征向量#这个特征向量本身也是一行数组#比如#矩阵A：# [[ 3 -1]# [-1  3]]#特征值a：[4. 2.]#特征向量b：[[ 0.70710678  0.70710678] [-0.70710678  0.70710678]]#array里存储的必须是同一类型的数据，list中的元素可以类型不同# Sorting  #信号子空间是由大的特征值对应的特征向量张成的子空间#噪声子空间是由小的特征值对应的特征向量张成的子空间#要找出噪声子空间，就要找小的特征值，要排序#接下来做的就是按照特征值的绝对值大小对特征向量排序eig_array = []for i in range(M):eig_array.append([np.abs(sigmai[i]),vi[:,i]])#先把特征值的绝对值和特征向量放到一行里，然后再把它们合并成一个元素放到一个list里#一共M(4)个接收机，协方差矩阵是4x4的，所以特征值和特征向量有4组#那么这个list有4个元素，每个元素的第一个子元素是特征值的绝对值eig_array[0]#第二个子元素是对应的特征向量#[[特征值1,[特征向量1[0],特征向量1[1],特征向量1[2],特征向量1[3]]],# [特征值2,[特征向量2[0],特征向量2[1],特征向量2[2],特征向量2[3]]],# [特征值3,[特征向量3[0],特征向量3[1],特征向量3[2],特征向量3[3]]],# [特征值4,[特征向量4[0],特征向量4[1],特征向量4[2],特征向量4[3]]]]eig_array = sorted(eig_array, key=lambda eig_array: eig_array[0], reverse=False)#sorted对所有可迭代对象进行排序#key是用来比较的元素，这个参数来自于可迭代的对象中，就是eig_array[0]#reverse=False，代表是从小到大排序的# Generate noise subspace matrixnoise_dimension = M - signal_dimension    #M是特征值总数，减去信号子空间的大小（1个信号发射机，信号子空间就对应1个特征值？）#就得到了噪声子空间的大小E = np.zeros((M,noise_dimension),dtype=complex)#E是4行3列的array，现在里面都是0#根据这个尺寸，应该是把3个噪声子空间对应的特征向量竖着存放，#每一列对应一个噪声子空间的特征向量for i in range(noise_dimension):     E[:,i] = eig_array[i][1] #i的范围是0，1，2，一开头是0#eig_array[i][1]，就是eig_array[0][1]，就是eig_array的第一个元素的第二个子元素，#经过排序，第一个元素就是特征值最小的那一对特征值+特征向量#第二个子元素就是这个特征向量了，所以E[:,0]里存的是特征向量#只不过E里是把特征向量竖过来存储的，所以是所有行（4行）和第一列是这个特征值向量#等所有3个噪声特征值对应的特征向量存完以后E就是噪声子空间了E = np.matrix(E)    theta_index=0for i in range(np.size(scanning_vectors, 1)):  #scanning_vectors里有4行361列，所以i的范围就是0～360           S_theta_ = scanning_vectors[:, i] #每行都要（对应4个接收机通道），列对应每个角度，这是阵列矩阵S_theta_  = np.matrix(S_theta_).getT() #从4行1列转为1行4列，暂时不知道为什么要转置？可能是为了下面的矩阵乘法？        ADORT[theta_index]=  1/np.abs(S_theta_.getH()*(E*E.getH())*S_theta_)#这个公式就是空间谱密度的公式theta_index += 1 #这个theta_index就跟i一样return ADORT

KerberosSDR代码笔记(5) 信号处理（采样时间延迟计算、相位差计算的2种方法、MUSIC算法）相关推荐

【MATLAB】数值计算：计算黄金分割比的N种方法（来自Matlab创始人Cleve Moler）
写作时间:2020-07-16 目录: 1.推荐一本书 2.开始第一课 3.总结一下正文 1.推荐一本书推荐一本书,<Numerical Computing with MATLAB>. ...
Shell脚本中计算字符串长度的5种方法
这篇文章主要介绍了Shell脚本中计算字符串长度的5种方法,来自于个人Shell脚本长期的开发经验,需要的朋友可以参考下有时在Linux操作系统中需要计算某个字符串的长度,通过查询资料整理了下目前S ...
Matlab 计算均方误差MSE的三种方法
Matlab 计算均方误差MSE的三种方法数据说明: ytest 测试集y,真实的y值,是一维数组: ytest_fit 基于测试集 x 预测的y值,是一维数组: test_error 是预测误差. ...
【MATLAB】数值计算：计算黄金分割比的N种方法
目录: 1.推荐一本书 2.开始第一课 3.总结一下正文 1.推荐一本书推荐一本书,<Numerical Computing with MATLAB>. 这本书是MATLAB 创始人 ...
matlab滤波器设计工具箱带阻滤波器,用matlab信号处理工具箱进行fir滤波器设计的三种方法...
用matlab信号处理工具箱进行fir滤波器设计的三种方法摘要介绍了利用 MATLAB 信号处理工具箱进行 FIR 滤波器设计的三种方法:程序设计法. FDATool 设计法和 SPTool 设 ...
KerberosSDR代码笔记(4) 接收机程序（滤波器、直流抑制、循环缓存）
接收机程序分为2部分,一部分在hydra_receiver.py中,另一部分在_receiver/C文件夹下,由几个c语言写的程序构成. 第二篇笔记中讲过: _receiver文件夹下还有个C文件夹, ...
Opencv3编程入门学习笔记（三）之访问图像像素的三种方法
访问图像像素的三种方法:指针访问,迭代器访问,动态地址访问.访问最快的为指针访问,以下算法在几毫秒,但指针访问容易造成内存泄漏:其次为迭代器访问:最后为动态地址访问. 以下程序是根据<OpenC ...
随机森林计算特征重要性_随机森林中计算特征重要性的3种方法
随机森林计算特征重要性 The feature importance describes which features are relevant. It can help with a better ...
java 文件checksum_计算文件Checksum的几种方法
回忆一下,自己是否在网站上下载文件时看到过Checksum这个东西,一串字符串? 比如,我们到Apache网站上去下载用于操作Excel的依赖包 - Apache POI,就可以看到checksum: ...

KerberosSDR代码笔记(5) 信号处理（采样时间延迟计算、相位差计算的2种方法、MUSIC算法）

KerberosSDR代码笔记(5) 信号处理（采样时间延迟计算、相位差计算的2种方法、MUSIC算法）相关推荐

最新文章

热门文章