スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

【MUSIC】SCILAB版 その2

SCILAB版を作成。

/////////////////////////////////////////////
//// MUSIC処理
//// 生成ファイル"music_data.dat"
//// scilab変数を保存したファイル
/////////////////////////////////////////////

clear

FILE_PASS = 'C:\scilab_works\MUSIC\';
FILE_NAME_R = 'array_data.dat';
FILE_NAME_W = 'music_data.dat';

load (FILE_PASS + FILE_NAME_R) // 受信信号 読込み

j = sqrt(-1); // jは複素数
pa = 3; // 推定する信号の個数

// 固有値(特異値)分解
[Uxx, Sxx, Vxx] = svd(Rxx);
sv = diag(Sxx);
Us = Uxx(:,1:pa);
Un = Uxx(:,pa+1:m);
// ステアリングベクトル(行列)の計算
th=[-90:1:90]';
A = exp(-j*2*%pi*([0:m-1]'*sin(th'*deg2rad))*d/wavelen)/sqrt(m);
// MUSICスペクトラムの計算
num = diag(A'*A);
En = Un'*A;
den = diag(En'*En);
doa = 10*log10(num./den);

// 推定結果の表示
// MUSIC Spectrum の表示
subplot(2,1,1);
plot(th, doa)
title('MUSIC DOA Estimator');
xlabel('Angle [deg.]')
ylabel('DOA Spectrum [dB]')
xgrid;

// 固有値の分布を表示
subplot(2,1,2);
plot(sv,'o')
title('Distribution of Eigenvalues');
xlabel('Index');
ylabel('Value')
xgrid;

save (FILE_PASS + FILE_NAME_W);
スポンサーサイト

【MUSIC】SCILAB版 その1

SCILAB版を作成。

/////////////////////////////////////////////
//// ANT受信信号生成
//// 生成ファイル"array_data.dat"
//// scilab変数を保存したファイル
/////////////////////////////////////////////

clear
//load sig_data // 変調波信号の読み込み

FILE_PASS = 'C:\scilab_works\MUSIC\';
FILE_NAME = 'array_data.dat';


///////////////// 信号生成 /////////////////
OS = 10; //オーバーサンプル数
n = 1:2048;
f1 = 1;
f2 = 2;
f3 = 3;
s1 = cos(2*%pi*f1/OS*n);
s2 = cos(2*%pi*f2/OS*n);
s3 = cos(2*%pi*f3/OS*n);
///////////////////////////////////////////


//
j = sqrt(-1); // jは複素数
m = 8; // m : アレー素子数
p = 3; // p:到来波数
fc = 1e9; // 搬送波の周波数(1GHz)
c = 3e8; // 伝搬速度
wavelen = c/fc; // 搬送波の波長
d = wavelen/2; // 素子間隔(半波長)
theta = [ 85; 40; -20]; // 入射角の設定
nn = 1024; // nn:データ数
sn1 = 5; // 信号1のSNR (dB)
sn2 = 5; // 信号2のSNR (dB)
sn3 = 5; // 信号3のSNR (dB)
SN = [ sn1; sn2; sn3 ]; // 各信号のSN比
sn = SN(1:p); // SNベクトルの定義
deg2rad = %pi/180; // ラジアンへの変換係数
//
// 信号源ベクトルの作成
tt = 1:nn;
//SS = [ s1(tt).'; s2(tt).'; s3(tt).'];
SS = [ s1(tt); s2(tt); s3(tt)];
S = SS(1:p,:);
// 観測雑音の生成
nr = rand(m,nn); //real
ni = rand(m,nn); //imag
U = nr+j*ni;
//
// SNの計算
//
Ps = S*S'/nn; // (ΣS)/nn
ps = diag(Ps); // 対角行列 ⇒ 各信号電力
refp= 2*10.^(sn/10); // SN真数変換 , "2"乗算は不明
tmp=sqrt(refp./ps); // SN調整して、各信号を電圧へ
S2 = diag(tmp)*S; // 対角化し、各SN値を各信号にのみ乗算
//S2 = [tmp(1,1)*S(1,:); tmp(2,1)*S(2,:); tmp(3,1)*S(3,:)]; と同等

//
// ステアリング行列の作成
//
a2=-j*2*%pi*([0:m-1]'*sin(theta'*deg2rad))*d/wavelen;
A = exp(a2); //行:ANT番号 列:信号種類(到来数)
//
// アレーアンテナ受信信号の計算
X = A*S2+U; //ステアリング*信号+ノイズ、信号内には各諸元をもつ波形の合成
//
// 相関行列と逆行列の計算
Rxx = X*X'/nn;
Rinv = inv(Rxx);
//
save (FILE_PASS + FILE_NAME, 'X', 'A', 'S2', 'U', 'Rxx', 'Rinv', 'm', 'nn', 'wavelen', 'd', 'deg2rad', 'theta', 'p');

直交復調(ゼロIF) シミュレーション

clear
clf()

//////////////// 共通パラメータ ///////////////////////////////////////////////////
n = 1:1024; //解析範囲
fs = 1000; //サンプリング周波数[MHz]
///////////////////////////////////////////////////////////////////////////////

Tini = 0;//%pi/4; //初期位相
f = 400; //RF周波数[MHz]
fl = 400; //LO周波数[MHz]


w = 2*%pi*f/fs;
wl = 2*%pi*fl/fs;

RF = cos(w*n + Tini);
LO_I = cos(wl*n);
LO_Q = cos(wl*n + %pi/2);

I = RF .* LO_I;
Q = RF .* LO_Q;

fig = figure(0);
subplot(121);
title("waveform");
xlabel("time");
ylabel("amplitude");
xgrid;
plot(n,RF)
plot(n,LO_I,'r')
plot(n,LO_Q,'g')
//zoom_rect(gca(),[0,0,20,1.5])

subplot(122);
title("waveform");
xlabel("time");
ylabel("amplitude");
xgrid;
plot(n,I)
plot(n,Q,'r')
//zoom_rect(gca(),[0,0,0.5,50])


///////////////////////////////////////////
// LPF
///////////////////////////////////////////
TAP_NUM = 64;

[wft,wff,fr]=wfir('lp',TAP_NUM,[0.1 0],'hn',[1 0]); //wft:インパルス応答、wff:周波数応答、fr:周波数グリッド

figure(1)
subplot(121)
title("周波数応答");
plot(fr,wff)
xgrid;

subplot(122)
title("インパルス応答");
plot(0:TAP_NUM-1,wft)
xgrid;


///////////////////////////////////////////
// LPF適用
///////////////////////////////////////////
I_LPF = convol(wft,I);
I_LPF_len = size(I_LPF,2); //畳込みの計算上、(信号長さ - 1) = (畳込み結果の長さ)
Q_LPF = convol(wft,Q);
Q_LPF_len = size(Q_LPF,2);

P = sqrt(I_LPF^2 + Q_LPF^2);

fig = figure(3);
subplot(121);
title("waveform");
xlabel("time");
ylabel("amplitude");
xgrid;
//plot(n,RF)
plot(0:I_LPF_len-1,I_LPF)
plot(0:Q_LPF_len-1,Q_LPF,'r')
//zoom_rect(gca(),[0,0,20,1.5])

subplot(122);
title("waveform");
xlabel("time");
ylabel("amplitude");
xgrid;
plot(0:I_LPF_len-1,P)
//plot(n,Q,'r')
//zoom_rect(gca(),[0,0,0.5,50])

途中まで

clear
delete("all")
clf()


TMP_D = "C\scilab_works\det_filter";
FILE_NAME = "test.txt";

//////////////// 共通パラメータ ///////////////////////////////////////////////////

FREQ_SAMP_ADC = 1000; //ADCサンプリングレート〔MHz〕


CAL_TIME = 20; //解析時間〔us〕◆ただし、FFTポイント数のため、2^Nへ変換している


TAP_NUM = 32; //フィルターのタップ数
///////////////////////////////////////////////////////////////////////////////


///////////////////////////////////////////
// LPF
///////////////////////////////////////////
[wft,wff,fr]=wfir('lp',TAP_NUM,[0.25 0],'hn',[1 0]); //wft:インパルス応答、wff:周波数応答、fr:周波数グリッド

figure(0)
subplot(121)
title("周波数応答");
plot(fr,wff)
xgrid;

subplot(122)
title("インパルス応答");
plot(0:TAP_NUM-1,wft)
xgrid;

///////////////////////////////////////////
// LPF ⇒ BPF化【入力信号のIQ分離】
// 「考え方①」実信号をIQ分離して、各々LPFへ入力すると考える。
// 「考え方②」実信号を複素BPF(IQそれぞれBPFをもつ)とする。
//
// 【メリット】
// 任意の実信号をLPF(FIR)へ入力した場合・・・
// LPFの動作周波数は、ADCサンプリング周波数と同じである必要がある。
// I成分BPF、Q成分BPFへ信号を2分配して入力した場合・・・
// 各BPF(FIR)の動作周波数は、ADCサンプリング周波数の半分でよい。
// ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
///////////////////////////////////////////

// e^(-jwn) = cos(wn) - j*sin(wn) , w = π/2とすれば、
// n = 0,1,2,3,4・・・のとき、1,-j,-1,j,1となり、IQが順番に出力されていることになる。
// ただし、"n"はサンプリング数なので、サンプリング周期(1/fs)で「1,-j,-1,j,1」は並んでいる
// IQごとに考えると、I「1,0,-1,0,1」,Q「0,-j,0,j,0」となる。
// つまり、IQを互いに"90度"ずらした"fs/2"で動作させ、最終的に加算すれば同じ結果となる。
// さらに、Iふたつ、Qふたつの計4つをすべて"90度"ずらせば、"fs/4"でよいことになる。

i_shift_def = [1 0 -1 0]; // 4サンプルで1周期とする
q_shift_def = [0 -1 0 1]; // 「i_shift_def」より、90度位相がずれている

I_shift = repmat(i_shift_def,1,TAP_NUM/4); //列方向に長さ拡張
Q_shift = repmat(q_shift_def,1,TAP_NUM/4);

bpf_I = wft .* I_shift; //LPFインパルス応答へシフト量だけ乗算
bpf_Q = wft .* Q_shift; //(時間軸での乗算 ⇒ 周波数軸での畳込み)

//FIR関数が出力するF特の周波数分解能と等しくするため
zeropad_num = size(wff,2)*2 - size(bpf_I,2); //FIR出力がナイキスト範囲なので2倍する
bpf_I = [bpf_I zeros(1,zeropad_num)];
bpf_Q = [bpf_Q zeros(1,zeropad_num)];

fft_bpf_I = abs(fft(bpf_I));
fft_bpf_Q = abs(fft(bpf_Q));

figure(1)
subplot(221)
title("周波数応答 BPF(I成分)");
plot((0:size(bpf_I,2)-1)/(size(bpf_I,2)-1),fft_bpf_I)
//plot((0:TAP_NUM-1)/(TAP_NUM-1),fft_bpf_I)
xgrid;

subplot(222)
title("周波数応答 BPF(Q成分)");
plot((0:size(bpf_Q,2)-1)/(size(bpf_Q,2)-1),fft_bpf_Q,'r')
//plot((0:TAP_NUM-1)/(TAP_NUM-1),fft_bpf_Q,'r')
xgrid;

subplot(223)
title("インパルス応答 BPF(I成分)");
plot(0:TAP_NUM-1,bpf_I(0:TAP_NUM-1))
xgrid;

subplot(224)
title("インパルス応答 BPF(Q成分)");
plot(0:TAP_NUM-1,bpf_Q(0:TAP_NUM-1),'r')
xgrid;

///////////////////////////////////////////
// 入力信号生成
///////////////////////////////////////////
fin = 50;//FREQ_SAMP_ADC * 0.5; // ? MHz
ts = 1/FREQ_SAMP_ADC;
N = 0:1:127;
freq_axis = (FREQ_SAMP_ADC/size(N,2))*N;

Ain = 2*%pi*fin*ts*N; // f*ts*N = (f/fs)*N
Sin = cos(Ain); //入力信号

fft_Sin = abs(fft(Sin));

figure(2)
subplot(121)
title("入力信号");
plot(N,Sin,'.-')
xgrid;

subplot(122)
title("FFT(入力信号)");
xlabel("Freq [MHz]");
ylabel("Power Spectrum")
plot(freq_axis,fft_Sin)
xgrid;

///////////////////////////////////////////
// IQ生成
//  入力信号をIQ分離するため
///////////////////////////////////////////

f = FREQ_SAMP_ADC * 0.125; //IQ周波数は、サンプル周波数 * 定数とする
ts = 1/FREQ_SAMP_ADC;
N = 0:1:127;
freq_axis = (FREQ_SAMP_ADC/size(N,2))*N;

A = 2*%pi*f*ts*N; // f*ts*N = (f/fs)*N
I = cos(A);
Q = sin(A); //(A-%pi*3/2);
amp0 = I + Q;
pow0 = sqrt(I^2 + Q^2);

fft_I = abs(fft(I));
fft_Q = abs(fft(Q));

//fft_RE = real(fft_Q);
//fft_IM = imag(fft_Q);

figure(3)
subplot(121)
title("IQ");
plot(N,I,'.-')
plot(N,Q,'.-r')
xgrid;

subplot(122)
title("FFT(IQ)");
xlabel("Freq [MHz]");
ylabel("Power Spectrum")
plot(freq_axis,fft_I)
plot(freq_axis,fft_Q,'r')
xgrid;


///////////////////////////////////////////
// 入力信号のIQ分離
///////////////////////////////////////////
S_I = Sin .* I; //I分離
S_Q = Sin .* Q; //Q分離
//POW_S_IQ = sqrt(S_I^2 + S_Q^2);


///////////////////////////////////////////
// フィルタ適用 (時間軸での畳込み計算)
///////////////////////////////////////////
conv_I = convol(wft,S_I);
conv_I_len = size(conv_I,2); //畳込みの計算上、(信号長さ - 1) = (畳込み結果の長さ)
conv_Q = convol(wft,S_Q);
conv_Q_len = size(conv_Q,2);

amp = conv_I + conv_Q;
pow = sqrt(conv_I^2 + conv_Q^2);

N2 = 0:size(conv_I,2)-1;
freq_axis_conv_IQ = (FREQ_SAMP_ADC/size(N2,2))*N2;
fft_conv_I = abs(fft(conv_I));
fft_conv_Q = abs(fft(conv_Q));

figure(4)
subplot(121)
title("フィルタと信号の畳込み");
xlabel("N sample");
plot(0:conv_I_len-1,conv_I,'.-')
plot(0:conv_Q_len-1,conv_Q,'.-r')
xgrid;

subplot(122)
title("FFT(convolution)");
xlabel("Freq [MHz]");
ylabel("Power Spectrum")
plot(freq_axis_conv_IQ,fft_conv_I)
plot(freq_axis_conv_IQ,fft_conv_Q,'r')
xgrid;


/////////////////////////////////////////////
//// フィルタ適用 (時間軸での畳込み計算)
/////////////////////////////////////////////
//conv_I = convol(wft,I);
//conv_I_len = size(conv_I,2); //畳込みの計算上、(信号長さ - 1) = (畳込み結果の長さ)
//conv_Q = convol(wft,Q);
//conv_Q_len = size(conv_Q,2);
//
//amp = conv_I + conv_Q;
//pow = sqrt(conv_I^2 + conv_Q^2);
//
//N2 = 0:size(conv_I,2)-1;
//freq_axis_conv_IQ = (FREQ_SAMP_ADC/size(N2,2))*N2;
//fft_conv_I = abs(fft(conv_I));
//fft_conv_Q = abs(fft(conv_Q));
//
//figure(4)
//subplot(121)
//title("フィルタと信号の畳込み");
//xlabel("N sample");
//plot(0:conv_I_len-1,conv_I,'.-')
//plot(0:conv_Q_len-1,conv_Q,'.-r')
//xgrid;
//
//subplot(122)
//title("FFT(convolution)");
//xlabel("Freq [MHz]");
//ylabel("Power Spectrum")
//plot(freq_axis_conv_IQ,fft_conv_I)
//plot(freq_axis_conv_IQ,fft_conv_Q,'r')
//xgrid;
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。