1. 自由研究に戻る
2. Menu表示

OctaveでMIMO、アレーアンテナの学習

　テキストのP64,【例8】　周波数選択性フェージングをの例を確認してみます。
シミュレーションのための仮定は
平均電力は遅延に対して指数関数的に減衰する仮定をします。
テキストには図が載ってないので、プロット図を書いてみます。

(1)遅延プロファイルの減衰量を10dBにして行ったシミュレーションです。

Fig.1は一様分布する位相を表わし、Fig.2はレイリー分布する振幅をしめします。振幅は複素振幅なので、plotはabs関数で大きさ表示しています。遅延100の時に第一到来波に対して10dB減衰する遅延プロファイルです。図のｘ軸delayは１delayあたり遅延量５に相当します(刻みが５ということ・・分かりずらい)。21delayで減衰量sqrt(10^(-1)/2)になるのは、信号振幅と平均電力の関係みたいです。
この位相と振幅でチャネル応答ベクトルを生成します。


Fig.3は送信信号です。BPSK変調し、ロールオフフィルタを通過させた信号です。 Fig.4はFig.3の信号に遅延を発生させ、先のチャネル応答ベクトルを掛算して作成した、フェージングを受けた信号になります。フェージングの影響を強く受けています。

(2)以下は、遅延プロファイルの減衰量を1000dBにして行ったシミュレーションです。

Fig.5はFig.1と同じ。Fig.6は遅延100の時に第一到来波に対して1000dB減衰する遅延プロファイルです。Fig.6のY軸がlogスケールになっていることに注意してください。Fig.2と減衰量が違いおおきな減衰です。delay21付近で減衰量sqrt(10^(-100)/2)になります。 ここのプロットは表示の関係でMatlabを使いました。Octaveではlog scaleのグラフで、y軸のスケールがうまく表示できませんでしたの、やむを得ず使いました。


Fig.7はFig.3と同じです。Fig.8はFig.7の信号に遅延を発生させ、チャネル応答ベクトルを掛算して作成した、フェージングを受けた信号になります。Fig.4に比べよりフェージングの影響が少なくなっています。しかしこの信号では受信できるかどうか不安です。

スペクトル拡散とRAKE受信

　テキストのP68のスペクトル拡散を確認してみます。
スペクトル拡散は周波数選択性フェージングの影響を軽減するために、先行波と遅延波を完全に分離する手法と言われています。完全に分離後、さらに先行波と遅延波をうまく合成できれば受信信号のSN比を改善できるというアイデアのようです。
また下記事項も重要と思われます。
・チャネル推定：SS（スペクトル拡散）信号を相関器に受信すると、遅延τのみあらず、通信路の位相・振幅を推定できる。
・RAKE合成：相関器によって分離分離されたパス毎のシンボルを合成することをRAKE合成といい、SN比を最大になるように合成するには最大比合成を用いる。

　 テキストの例題を見てみます。
プログラム例10　SS信号の生成、送信と受信
（未実施）
プログラム例11　RAKE合成のシミュレーション
（未実施）

OFDM

　　 受信側に複数のアンテナを設置して、安定した受信を実現する。受信側で行う空間ダイバーシティには３通り考えられる。
　 （１）選択ダイバーシティ　
　 複数のアンテナのうち、どれかを選択して使う。
　 （２）等利得合成　
　 複数のアンテナの位相を合わせて合成する。アンテナの信号強度が落ち込んでいるとSNが劣化する。
　 （３）最大比合成　
　 位相と振幅を合わせて合成する。SNが最大になるように合成する　

　送受信モデルを考える。送受信アンテナ間の応答をここではチャンネル応答と呼んでいる。すなわち送受信の伝達関数（チャネル応答）を求め、伝達関数のパラメータを同定する問題に帰着する。

重み係数W（は「複素共役して掛けるもの」と考えておくと後々理解しやすい）
　出力 Y[k] = W^TX[k]
　W=[w1 ; w2 ; ;wn]
　最大比合成された信号は、各アンテナのSN比の合計に等しい

・最大比合成重み係数の基本的考え方は２通り
　チャネル推定（チャネル応答を用いる）
　　チャネル応答ベクトルｈはそのまま最大比合成重み係数になる
　　ｈの共役複素を受信データに掛けるとチャネルの位相回転はキャンセルされる
　固有ベクトルを用いる方法
　

最大比合成の確認

　テキストから少しずれますが、参考文献（９）を元に、SFBC(space-frequency block coding)という送信ダイバーシティを眺めてみます。２アンテナポートの場合、LTEではSFBCダイバーシティが用いられます。最初のアンテナポートには連続する変調シンボルd(i)とd(i+1)とが互いに隣接するサブキャリアにマッピングされます。第2のアンテナポートには-conj (d(i+1))とconj (d(i))とが同一（最初のアンテナポートと同じという意味）の隣接サブキャリアで送信されます。ただしconj ()は複素共役を表わします。
　3GPPでは TS36.211 6.3.3.3及び6.3.4.3でSFBCが定義されて、次のように表わしています。

　3GPP TS36.211 6.3.3.3でのlayer mapping
　x⁽⁰⁾(i)=d⁽⁰⁾(2i)
　x⁽¹⁾(i)=d⁽⁰⁾(2i+1)
　
　3GPP TS36.211 6.3.4.3でのprecoding
　y⁽⁰⁾(2i)=(1/sqrt(2)){Re(x⁽⁰⁾(i) + j Im(x⁽⁰⁾(i)} = S1
　y⁽¹⁾(2i)=(1/sqrt(2)){-Re(x⁽¹⁾(i) + j Im(x⁽¹⁾(i)} = -S2^*
　y⁽⁰⁾(2i+1)=(1/sqrt(2)){Re(x⁽¹⁾(i) + j Im(x⁽¹⁾(i)} = S2
　y⁽¹⁾(2i+1)=(1/sqrt(2)){Re(x⁽⁰⁾(i) - j Im(x⁽⁰⁾(i)} = S1^*

参考文献（９）P36より

サブキャリア	アンテナ１	アンテナ２
2k	S1	-S2*
2k+1	S2	S1*

この後は、MRCを用いて受信するようです。

　 空間ダイバーシティは、アンテナを空間的に十分離して設置していました。
アダプティプアレーは、反対にアンテナ間を波長間隔（半波長も含めて）に設置、受信信号に相関をもたせています。アンテナの指向性も構成できます。
　 受信信号のモデルには下記がポイント
　　 平面波の仮定
　　 経路長差
アンテナ間隔をd、角度θで信号が到来する。ちょうど光学部品の回折格子のような振る舞いをします。
経路長差はdsinθ、位相差はφ= 2π/λ・dsinθ

ステアリングベクトルの導出
　 steering vector（steerは操縦する、案内するという意味）
X[k]=A(θ)S[k]　： X[k]は受信信号ベクトル、A(θ)はステアリングベクトル、S[k]は信号

相関行列
アンテナ素子間の相関特性（コヒーレンス）を行列で表わす。(k,l)成分が素子kと素子lの間の受信信号の相関を表わすようにした行列である。
Rxx=E[X(t)・X(t)^H]
E[・]は期待値を意味します。期待値は平均と類義語です。


ウィナーの解
　 希望信号を取り込みながら、干渉信号を打ち消す。このような重み係数を求めるにはどうすのだろうか。このような状態を完全に解くことは出来ないが、最もよい解をえることは出来る。この解がウィナー解である。
それは評価関数を最小にする、最適制御と同じ考えである。２乗平均誤差を最小にする、２乗誤差最小化（MMSE）基準を満足するような解と同じである。

%ウィーナー解の計算
相関行列Rxx
参照信号ベクトルｒ
アンテナ入力信号X
ウィーナー解 w=inv(Rxx)*X*conj(r)
　　 X*conj(r)は参照信号と入力信号ベクトルとの間の相関ベクトル
アンテナ出力は y=w'*X

MMSEアダプティブアンテナ
MMSE基準（または規範）のアダプティブアンテナは、受信側で用意する参照信号と実際のアレー出力徒の差を最小にすることによって最適なウェイトを決定する。アダプティブヌルステアリングと同時にアダプティブビームフォーミングを行うことが出来るが、参照信号が必要になる。
・LMS(Least Mean Square)法 --ウィナー解の最適化アルゴリズムの一つ、最急降下法
・SMI(Sample Matrix Inversion)法-- ウィナー解の理論式を直接解く直接解法
・RLS(Recursiv Least Square)法-- SMIのようにサンプル値を用いて相関値を推定、かつウェイトを逐次更新

アンテナ指向性の描画
　 アンテナ数、波長で正規化したアンテナ間距離、および重み係数が分かれば、アンテナの指向性を調べる事が出来る。

出力SN比の計算法
　 アダプティブアレーアンテナによって出力SN比は入力SN比のアンテナ数倍大きくなる。ただし干渉波を打ち消すことを同時に行っている場合はそうならない。

逐次的な重み係数うの推定法
LMS法
RLS法
重み係数のブラインド推定 CMA法

ウィナー解に基づくアダプティブアレーアンテナ

　テキストのP122,【例14】　ウィナー解に基づくアダプティブアレーアンテナの確認をして見ます
　ここでの方法は、直接解法なのでSMI法に当てはまるようです。
　設定項目は
　　 　希望信号到来角度=-20
　　 　干渉信号到来角度=30
　　 　アンテナ数=4
　　 　アンテナ間隔=0.5
　　 　SN

　希望到来信号 S(1,:)
　干渉到来信号 S(2,:)
　参照信号をr=S(1,:)として、wienerの解で重み係数を求めたあとで、その重み係数を用いてデータ信号を受信する。
　 　 （１）送信データ（BPSK）および、アンテナ入力信号生成

Fig.1 送信データです（BPSK）
Fig.2はステアリングベクトルと複素ガウス雑音を加え、アンテナ入力信号Xを計算する。4アンテナ分をプロット。

　 （２）ウィーナーの解に基づくアンテナ出力

Fig.1はウィーナーの解に基づき重み係数を求め、アンテナ出力信号を生成しています。
Fig.2は元データです。
%ウィーナー解の計算
相関行列Rxx
参照信号ベクトルｒ
アンテナ入力信号X
ウィーナー解 w=inv(Rxx)*X*conj(r)
　 アンテナ出力は y=w'*X
入力SNは10dBで出力SNは15.58dBであった。理論的にはアンテナ４本なので6dB高くなる。
それでは
アンテナパターンを見てみます。テキストでは横軸を角度にしていますが、極座標で見てみます。
ただし、アンテナは等間隔で一列に並んでいる場合ですので、並べ方を変えれば違うパターンになるでしょうから注意してください。

Octaveには極座標でのplot関数polarが用意されています
例えば、polar(doa*pi/180,10*log10(P_theta)+50)のように使います

（２−１）アンテナ数４の場合（【例14】本例題の場合）

Fig3はアンテナ数４の場合です。
希望信号到来角度=-20(赤色ライン)、干渉信号到来角度=30（緑色ライン）
信号強度は相対値です。-40,-20とか負号がついていますが、Octaveが付けているもので、方向（角度）を意識して描いているのでしょう。信号強度が負なのではありません。
図において水平方向（X軸）の右側は０度、左側が−180度になりますので、４本のアンテナは垂直軸（Ｙ軸）方向に並んでいることになります。4本の各アンテナは、それぞれ図の右半面と左半面対象にアンテナパターンが描かれています。

（２−２）アンテナ数２の場合（【例14】のアンテナ数を２にする）

Fig4はアンテナ数２の場合です。

（２−３）アンテナ数８の場合（【例14】のアンテナ数を８にする）

Fig5はアンテナ数８の場合です。
希望信号到来角度=-20度(赤色ライン)、干渉信号到来角度=30度（緑色ライン）
何度かシミュレーションを実行してみますと、うまくいかない（希望波をとらえられない）場合もあるようです。 アンテナを増やせばうまくゆくとも限らないような気がします。

LMS法/RLS法

逐次的な重み係数の推定法としてMMSE基準によって最適解を求めるLMS法とRLS法が紹介されています。 LMS法はサンプルごとに重み係数を更新していきますが、LMS更新式の収束を左右するステップサイズの取り方が難しい。ステップサイズが大きいと収束は早いが不安定になり、小さいと収束は遅いが安定となる。
ここで考えられたのがRLS法で、ステップサイズを適応的に変えていく方法です。

LMS法のアンテナパターン

（３）テキストのプログラム【例15】P124のアンテナパターンを描いてみます。アンテナ数２，４として希望波、干渉波の到来角度を変化させてみてみます。


Fig.lms1 アンテナ数２
希望信号到来角度=-20度(赤色ライン)
干渉信号到来角度=30度（緑色ライン）

Fig.lms2 アンテナ数２
希望信号到来角度=0度(赤色ライン)
干渉信号到来角度=30度（緑色ライン）


Fig.lms3 アンテナ数２
希望信号到来角度=-20度(赤色ライン)
干渉信号到来角度=0度（緑色ライン）

Fig.lms4 アンテナ数４
希望信号到来角度=-20度(赤色ライン)
干渉信号到来角度=30度（緑色ライン）


Fig.lms5 アンテナ数４
希望信号到来角度=0度(赤色ライン)
干渉信号到来角度=30度（緑色ライン）

Fig.lms6 アンテナ数４
希望信号到来角度=-20度(赤色ライン)
干渉信号到来角度=0度（緑色ライン）

LMS法のアンテナパターン、干渉信号が２波の場合

プログラム【例14１】〜【例16】は干渉波を増やすことが出来ます。
干渉波を2つにしてみます。

Fig.lms7 アンテナ数４
希望信号到来角度=-20度(赤色ライン)
干渉信号到来角度=0度（緑色ライン）
干渉信号到来角度=30度（黄色ライン）

Fig.lms8 アンテナ数８
希望信号到来角度=-20度(赤色ライン)
干渉信号到来角度=0度（緑色ライン）

RLS法のアンテナパターン

（４）テキストのプログラム【例16】P125のアンテナパターンです。
wiener法、LMS法、RLS法は重み係数を求める方法が違いますが、このシュミレーション条件ではアンテナパターンに大きな違いが見えないようです。
LMS法の場合と同じように希望波の到来角度を変化させたときのアンテナパターンをRLS法の場合に見てみます。先ほどと同じように角度によって”死角”（不得意角度）のようなものがあるようです。（残念ですがこれはアンテナを一列に配置した場合ですのでそれ以外の配置に当てはまるかどうかは議論されていません）


Fig.RLS1 アンテナ数４
希望信号到来角度=-20度(赤色ライン)
干渉信号到来角度=30度（緑色ライン）

Fig.RLS2 アンテナ数４
希望信号到来角度=-40度(赤色ライン)
干渉信号到来角度=30度（緑色ライン）


Fig.RLS3 アンテナ数４
希望信号到来角度=-50度(赤色ライン)
干渉信号到来角度=30度（緑色ライン）

Fig.RLS4 アンテナ数４
希望信号到来角度=-60度(赤色ライン)
干渉信号到来角度=30度（緑色ライン）

CMA法を用いてみます

CMA法によるアダプティブアレーアンテナは、参照信号を不要とする重み係数推定法です。
テキストの例17です。
（５）CMA(Constant modulus algorithm)法によるアンテナ出力

Fig.3CMA法によるアンテナ出力です。
Fig.4は元データです。SNは14.43とwiener法に比べ低い値です。
CMA法は、アダプティブアレーアンテナ出力の包絡線が一定となることを制御基準とする。PSK変調には適用できると言いますが、OFDMには適用できなくなるということでしょう。

　 送信側にも複数のアンテナを用いることで、マルチパスフェージング環境で通信路容量を増やせる。MIMOも空間ダイバーシティの一つ。
X=H*U+雑音
X=MIMOの受信信号
H=チャネル応答行列
U=送信信号行列

H=MYcompNoise([N_R,N_T],1);%チャネル応答ベクトル 平均電力１の複素正規乱数
[V_R,LambdaM,V_T]=svd(H);%特異値分解

１、特異値分解に基づくマルチモード伝送の例です
テキスト　プログラム【例18】の確認です。SNは20dBで計算。
特異値分解では U=VtS：Sは信号行列
アンテナ N_R=2,N_T=2の 2×2MIMO, Ndata=1000,BPSKデータです。
（１）信号系列を作る。ランダムなBPSK信号です。

Fig.1はアンテナ１の信号系列です。ランダムデータのBPSK信号です（100シンボルまで表示）
Fig.2はアンテナ２の信号系列です。

（２）送信信号系列を作る。U=V_T*S

Fig.3,4は信号系列の実部と虚部です。青色：アンテナ１、緑色：アンテナ２です。実部を見ると、アンテナ１の値が平均して大きい。虚部はアンテナ２が大きいように見えます。
　 　
（３）受信信号系列Yを求める

Fig.5,6はアンテナ１と２の受信信号系列です。アンテナ１の信号が2倍の振幅で受信されています。これはMIMOの性質なのでしょうか。テキストには、それぞれの伝送路のSN比は異なるという記述があります。
Y=V_R'*X; %ただし[V_R,LambdaM,V_T]=svd(H);
（４）受信信号系列を復調する。

Fig.7,8は復調したデータです。Fig.1,2（下の図）と比較してみます。とりあえず目視比較で、一致していそうです。


次に

２、OFDM-64QAM信号で確認してみます
雰囲気を味わうため、計算を簡単にします。遅延なし、フェージングなし、GIなしで考えてみます。
チャンネル推定様にM系列を先頭に付加します。OFDMの項で確認した内容と同じです。SN40dBを設定。
（遅延とフェージングが無いと位相が回らないようなので、後ほど入れて検討してみます）

（１）M系列を含んだ信号系列を作る　アンテナ１

Fig.9-1は信号の先頭につけた”チャネル推定用系列”のM系列信号です。
Fig.10-1は変調信号全体ですが、65408サンプルまではチャネル推定用系列のM系列が付加されています（長さは511×128=65408）、それ以降はデータです。正確な。データ長は128×1000=128000ですので　全体で193408です。

（２）M系列を含んだ信号系列を作る　アンテナ２


（３）送信信号系列Uを作る　アンテナ１と２

Fig.11-1,11-2は送信信号系列はU=V_T*Sですが、アンテナ１とアンテナ２ではM系列に相当する部分の振幅が異なっています。特異値分解されたV_Tの性質なんでしょうか。？これは特徴のようです。

（４）MIMO受信信号Xを求める　アンテナ１と２

Fig.12-1, 12-2はMIMO受信信号X=H*V_T*S+複素ガウスノイズです。

（５）受信信号行列Yを求める　アンテナ１　

Fig.13-1,14-1はY=V_R'*X;%受信信号出力（実部と虚部）です。受信信号Xに受信重みV_Rを掛けてます。
これで特異値分解によるMIMO分離が出来ているようです。

コンステレーションを見てみます。FFTしてコンステレーションを描きます。
scatterplot(FFTout1(1,512:1511),1,0,'b.');シンボルは512〜1511に有ります。
振幅が大きくなってます、振幅補正が出来てない為です。後で振幅補正します。


（６）受信信号行列Yを求める　アンテナ２


（７）OFDMのためのチャネル推定を行う　アンテナ１、アンテナ２
チャネル推定後、位相とゲインの補正が必要です。
%受信機　チャネル推定　MIMOの場合 chEstSeq_rx1=FFTout1(:,1:seqLength);
phaseEst1=chEstSeq_rx1*mseq;
phaseEst12=phaseEst1;
phaseEst1=phaseEst1./abs(phaseEst1);%チャネル推定値

%比較のためもとの振幅１、位相遅れ無を確認する
phaseOrg1=chEstSeq1*mseq;
phaseOrg2=chEstSeq2*mseq;
%受信シンボルの位相補償とQAMの復調 振幅も補償する
%antenna1
phaseCompMat1=phaseEst1*ones(1,NofdmSymbol); %位相補償の値
rxSymbol12=FFTout1(:,(seqLength+1:end)).*conj(phaseCompMat1);
phaseEst13=phaseEst12./phaseOrg1;
gainCompMat1=phaseEst13*ones(1,NofdmSymbol); %振幅補償の値
rxSymbol13=rxSymbol12./abs(gainCompMat1);ここで振幅補償している。
rxSymbol1=round(rxSymbol13);
rxData1=qamdemod(MYvec(rxSymbol1),M);

Fig15-1は位相補償の図です、位相遅れはないようです。振幅は１で正規化しています。
Fig15-2は振幅補償の図です。振幅は2倍近くになっていて、位相遅れがないことがわかります。


振幅補正後のコンステレーションを描いてみます。予定道理の６４QAMコンステレーションです。
scatterplot(rxSymbol13(1,:),1,0,'b.');


（８）受信シンボル　アンテナ１


（９）受信シンボル　アンテナ２


（１０）受信（復調後）データ、送信元データ比較　アンテナ１


（１１）受信（復調後）データ、送信元データ比較　アンテナ２

３、OFDM-64QAM信号 フェージングありの場合を確認してみます
遅延、フェージングを考慮し、GIを含めて計算してみます。
チャンネル推定様にM系列を先頭に付加します。OFDMの項で確認した内容と同じです。SN50dBを設定。
N_R=2;%受信アンテナ数
N_T=2;%送信アンテナ数
Nsub=128;% サブキャリア数
NofdmSymbol=1000;%OFDMシンボル数
Ns=Nsub*NofdmSymbol; %データシンボル数
N_GI=16; %GIの長さ
%フェージング設定
delay=[0;5];% 遅延τの設定
Tau1=30;
Mitigation_db=-10;
channelResponse1=MYdelayProfile(delay,Tau1,Mitigation_db);
channelResponse2=channelResponse1;

（１）M系列を含んだ信号系列を作る。ＧＩも付加します　アンテナ１
%送信機；データ信号生成 OFDM　64QAM
M=64;
data1=round(rand(Ns,1)*(M-1));%antenna1
qamSymbol1=qammod(data1,M);
data2=round(rand(Ns,1)*(M-1));%antenna2
qamSymbol2=qammod(data2,M);
SPCout1=reshape(qamSymbol1,Nsub,NofdmSymbol);%　S/P変換 パラレルに変換 antenna1
SPCout2=reshape(qamSymbol2,Nsub,NofdmSymbol);%　S/P変換 パラレルに変換 antenna2
%送信機　チャネル推定用系列
mseq=MYmseq([0;0;0;1;0;0;0;0;1]); % M系列の発生
seqLength=length(mseq); %チャネル推定系列長
chEstSeq1=ones(Nsub,1)*mseq.';%antenna1
chEstSeq2=chEstSeq1;%antenna2は１に同じ
%送信機 IFFT　（GI付加)
IFFTout1=ifft([chEstSeq1,SPCout1]);%アンテナ１
GI1=IFFTout1((end-N_GI+1):end,:);
IFFTout2=ifft([chEstSeq2,SPCout2]);%アンテナ２
GI2=IFFTout2((end-N_GI+1):end,:);
sOFDM1=MYvec([GI1;IFFTout1]);
sOFDM2=MYvec([GI2;IFFTout2]);
sOFDM =[sOFDM1,sOFDM2];%antenna1と2を束ねる

Fig.20-1は信号の先頭につけた”チャネル推定用系列”のM系列信号です。
Fig.21-1は変調信号全体ですが、65408サンプルまではチャネル推定用系列のM系列が付加されています（長さは511×128=65408）、それ以降はデータです。正確なデータ長は128×1000=128000とＧＩ長128×16=2048ですので　全体で193408＋2048=195456です。

（２）M系列を含んだ信号系列を作る。ＧＩも付加します　アンテナ２


（３）送信信号系列Uを作る　アンテナ１と２
周波数選択性フェージングを考慮してみます。アンテナ１，２とも同じフェージング特性としています。
H=MYcompNoise([N_R,N_T],1);%チャネル応答ベクトル 平均電力１の複素正規乱数
[V_R,LambdaM,V_T]=svd(H);%特異値分解
U=V_T*sOFDM.';
U1=U(1,:);
U2=U(2,:);
UU1=MYfSelFading(U1,delay,channelResponse1);
UU2=MYfSelFading(U2,delay,channelResponse2);
UU=[UU1.';UU2.'];
SNRdB=50;%
Pn=10^(-SNRdB/10);
X=H*UU+MYcompNoise(size(UU),Pn);%アンテナ入力信号
Y=V_R'*X;%受信信号出力

Fig.22-1,22-2は送信信号系列はU=V_T*Sですが、アンテナ１とアンテナ２ではM系列に相当する部分の振幅が異なっています。特異値分解されたV_Tの性質なんでしょうか。？これは特徴のようです。

（４）MIMO受信信号Xを求める　アンテナ１と２

Fig.23-1, 23-2はMIMO受信信号X=H*V_T*S+複素ガウスノイズです。

（５）受信信号行列Yを求める　アンテナ１　


（６）受信信号行列Yを求める　アンテナ２


（７）OFDMのためのチャネル推定を行う　アンテナ１、アンテナ２
%ここは　OFDM受信機
%rxSPCOut1=reshape(Y(1,:),Nsub,NofdmSymbol+seqLength); %S/P
rxSPCOut1=reshape(Y(1,:),Nsub+N_GI,NofdmSymbol+seqLength); %S/P
rxSPCOut1(1:N_GI,:)=[];%GI除去
FFTout1=fft(rxSPCOut1);
%rxSPCOut2=reshape(Y(2,:),Nsub,NofdmSymbol+seqLength); %S/P
rxSPCOut2=reshape(Y(2,:),Nsub+N_GI,NofdmSymbol+seqLength); %S/P
rxSPCOut2(1:N_GI,:)=[];%GI除去
FFTout2=fft(rxSPCOut2);
%受信機　チャネル推定　MIMOの場合
chEstSeq_rx1=FFTout1(:,1:seqLength);
phaseEst1=chEstSeq_rx1*mseq;
phaseEst12=phaseEst1;
phaseEst1=phaseEst1./abs(phaseEst1);%チャネル推定値
chEstSeq_rx2=FFTout2(:,1:seqLength);
phaseEst2=chEstSeq_rx2*mseq;
phaseEst22=phaseEst2;
phaseEst2=phaseEst2./abs(phaseEst2);%チャネル推定値

振幅、および位相の変化が発生しています。

（８）受信シンボル　アンテナ１
%受信シンボルの位相補償とQAMの復調 振幅も補償する
%antenna1
phaseCompMat1=phaseEst1*ones(1,NofdmSymbol); %位相補償
rxSymbol12=FFTout1(:,(seqLength+1:end)).*conj(phaseCompMat1);
phaseEst13=phaseEst12./phaseOrg1;
gainCompMat1=phaseEst13*ones(1,NofdmSymbol); %振幅補償
rxSymbol13=rxSymbol12./abs(gainCompMat1);
rxSymbol1=round(rxSymbol13);
rxData1=qamdemod(MYvec(rxSymbol1),M);


（９）受信シンボル　アンテナ２


（１０）受信（復調後）データ、送信元データ比較　アンテナ１

うまく受信されている様子です。

（１１）受信（復調後）データ、送信元データ比較　アンテナ２

　Zero ForcingによるMIMOを考えます。チャネル応答にその逆システム（逆行列）を掛けてチャネル応答を打ち消す方法は、ノイズ増幅の問題はあるものの、もっとも手っ取り早いMIMOの手法です。
X=H*U+雑音
X=MIMOの受信信号
H=チャネル応答行列
U=送信信号行列

　Zero ForcingによるMIMOを考えます。チャネル応答にその逆システム（逆行列）を掛けてチャネル応答を打ち消す方法は、ノイズ増幅の問題はあるものの、もっとも手っ取り早いMIMOの手法です。
X=H*U+雑音
X=MIMOの受信信号
H=チャネル応答行列
U=送信信号行列

Nsub=20;% サブキャリア数（行に割り当て）
NofdmSymbol=1000;%変調シンボル数（列に割り当て）
Ns=Nsub*NofdmSymbol; %64QAMシンボル数
N_GI=5;% GI(gard interval) length

%送信機；データ信号生成 64QAM
M=64;
data=round(rand(Ns,1)*(M-1));
qamSymbol=qammod(data,M);
SPCout=reshape(qamSymbol,Nsub,NofdmSymbol);%　S/P変換

y=fft(SPCout);%SPCoutが配列なので列方向にfftする

%送信機　チャネル推定用系列
mseq=MYmseq([0;0;0;1;0;0;0;0;1]); % M系列の発生
seqLength=length(mseq); %チャネル推定系列長 511
chEstSeq=ones(Nsub,1)*mseq.';

%送信機 IFFTおよびGIの付加
IFFTout=ifft([chEstSeq,y]); %チャネル推定は先頭の列に付ける

GI=IFFTout((end-N_GI+1):end,:);
sSCFDMA=MYvec([GI;IFFTout]);%GIを先頭の行に付ける
figure(1)
subplot(2,1,1)
plot([real(sSCFDMA(1:(N_GI+Nsub)*511)),imag(sSCFDMA(1:(N_GI+Nsub)*511))])

grid 'on'
title('Fig.1 sSCFDMA GI+chEst ,Blue=real,Green=imag') %chEst部のGIも含む

subplot(2,1,2)
plot([real(sSCFDMA(1:(Nsub*(NofdmSymbol+511)+N_GI*511))),imag(sSCFDMA(1:(Nsub*(NofdmSymbol+511)+N_GI*511)))])
title('Fig.2 (GI +chEst +  data) Blue=real,Green=imag') 

%フェージング設定
delay=[0;5];% 遅延τの設定
Tau1=30;
Mitigation_db=-10;
channelResponse=MYdelayProfile(delay,Tau1,Mitigation_db);

%周波数選択制フェージング通信路
chOut=MYfSelFading(sSCFDMA,delay,channelResponse);
SNRdb=100; %SNRの設定
Pn=10^(-SNRdb/10)/Nsub;
chOut=chOut+MYcompNoise(size(sSCFDMA),Pn);%雑音の付加

figure(2)
subplot(2,1,1)
plot([real(chOut(1:1000)),imag(chOut(1:1000))])
grid 'on'
title('Fig.3 [ real(chOut(1:1000)),imag(chOut(1:1000)) ]') 
 
subplot(2,1,2)
plot([real(chOut(1:(Nsub*(NofdmSymbol+511)+N_GI*511))),imag(chOut(1:(Nsub*(NofdmSymbol+511)+N_GI*511)))])
title('Fig.4 (ChEst + GI + data + noise)Blue=real,Green=imag') 

%受信機
rxSPCOut=reshape(chOut,Nsub+N_GI,NofdmSymbol+seqLength); %S/P
rxSPCOut(1:N_GI,:)=[];%GI除去
FFTout=fft(rxSPCOut);

%受信機　チャネル推定
chEstSeq_rx=FFTout(:,1:seqLength);
phaseEst=chEstSeq_rx*mseq;
phaseEst2=phaseEst;
phaseEst=phaseEst./abs(phaseEst);%チャネル推定値

%比較のためもとの振幅１、位相遅れ無を確認する
phaseOrg=chEstSeq*mseq;
%phaseOrg=phaseOrg./abs(phaseOrg);%チャネル推定値
subplot(2,1,2)
plot([real(phaseOrg),imag(phaseOrg)])
%ylim([-2 2])
title('Fig.6 [real(originalPhase),imag(originalPhase)]') 

%受信シンボルの位相補償とQAMの復調 振幅も補償する
phaseCompMat=phaseEst*ones(1,NofdmSymbol); %位相補償
rxSymbol2=FFTout(:,(seqLength+1:end)).*conj(phaseCompMat);

phaseEst3=phaseEst2./phaseOrg;
gainCompMat=phaseEst3*ones(1,NofdmSymbol); %振幅補償

rxSymbol3=rxSymbol2./abs(gainCompMat);
rxSymbol=round(rxSymbol3);

rxData=qamdemod(MYvec(rxSymbol),M);

rxDataIFFT=ifft(rxData);
rxSymbolIFFT=ifft(rxSymbol);