نقدر إنتاجية قناة MIMO (يتم تضمين خوارزمية صب الماء)


مقدمة


في صيف عام 2016 ، من حدث معروف ، تمكن خادمك المتواضع ، من بين مجموعة من الطلاب الآخرين ، من حضور محاضرات البروفيسور مارتن هارارت حول ميمو ، والتي عقدها في إطار برنامج الماجستير الدولي "الاتصالات ومعالجة الإشارات" . لكن لسوء الحظ ، مرضت أسبوعًا ونصفًا من الأسبوعين - ومن ثم فقد خرج عدد من الموضوعات عن نطاق تفهمي ... ومع ذلك ، بعد مرور بعض الوقت ، أصبح تحليل أساسيات MIMO هوايتي - وليس ترك القضية غير مكتملة.


شيئًا فشيئًا ، تطور كل هذا ليصبح سلسلة من الحلقات الدراسية المجردة الصغيرة ، والتي ربما تكون من الخطأ عدم مشاركتها. واليوم ، تكريما ليوم الاتصالات ، أود أن أناقش معك موضوع عرض النطاق الترددي لقناة MIMO - موضوع بسيط ، ولكن لا يزال يسبب بعض الصعوبات للطلاب (وليس الطلاب فقط).


قد يبدو للناس غير المعنيين أن زيادة عدد هوائيات الاستقبال والإرسال في إطار التقنية المذكورة أعلاه تزيد من عرض نطاق النظام بنفس المقدار: على سبيل المثال ، إذا وضعت هوائيات على جانب الاستقبال و 2 هوائيات على جانب الإرسال (MIMO 2x2) 2 مرات. ولكن هل هذا صحيح حتى من الناحية النظرية؟ دعنا نحاول معرفة ذلك!


يمكن العثور على نسخة أكثر رسمية باللغة الإنجليزية على الرابط وفي مستودع جيثب الخاص بي.
في هذه المقالة ، لن نفكر في مسائل ارتباط الهوائي ومشاكل التنفيذ الأخرى. نحن نحصر أنفسنا في نظرية مقطرة - بالنسبة للمبتدئين.

تلقى نموذج إشارة



قبل أن نبدأ الحديث عن النطاق الترددي ، سنتعامل أولاً مع الوصف الرياضي للإشارة المستلمة. يجب أن يؤخذ هذا الجزء بعناية فائقة ، لأن الكثير سيأتي من هذه الصيغة. لذلك يمكن تمثيل الإشارة المستلمة على النحو التالي:


\ mathbf {y} = \ sqrt {\ frac {P} {M_T}} \ mathbf {H} \ mathbf {s} + \ mathbf {n} \ qquad (1)

حيث P - قدرة الارسال ، M_T - عدد هوائيات الإرسال ، \ mathbf {s} - الشخصيات المنقولة \ mathbf {n} - الضوضاء المضافة ، و \ mathbf {H} - مصفوفة معاملات الإرسال للقناة (في الواقع ، عملية الخبو).


يمكن أيضًا إرسال الإشارة المرسلة بمزيد من التفصيل:


s_i = \ gamma_i d_i \ quad i = 1،2 ، .. M_T \ qquad (2)

حيث دي - واحدة من إشارات المعلومات ( E \ {\ mathbf {d} \ mathbf {d} ^ H \} = M_T ) ، و \ gamma_i - تضخيم مسير معين من انتشار موجة EM (كسب المسير).


دعونا نذكر ما هو تعزيز المسار:
إن تقوية مسار الانتشار (أو وزن الهوائي) يعني توزيع قدرة الخرج بما يتناسب مع "شدة" مسار معين. بعبارة أخرى ، نريد تخصيص المزيد من القدرة للقنوات الجيدة (مسارات الانتشار) وخفض الطاقة للقنوات السيئة.


أوزان الهوائي محدودة بعدد هوائيات الإرسال:


\ sum ^ r_ {i = 1} \ gamma_i = M_T \ qquad (3)

حيث ص هي رتبة مصفوفة القناة.


الحديث عن هذا الأخير.


مصفوفة البعد \ mathbf {H} يعوض M_R \ الأوقات M_T حيث M_R - عدد هوائيات الاستقبال.
لعدة قياسات زمنية ، ستبدو القناة كما يلي:



للرجوع اليها:
ربما بالنسبة للحسابات والنماذج الأكثر تعقيدًا ، سترغب في استخدام واحدة من أكثر الأدوات شعبية لذلك - MatLab . في هذه الحالة ، تجدر الإشارة إلى أنه يتم استخدام بنية بيانات مختلفة قليلاً هناك: الصفوف هي قياسات زمنية (لقطات) ، وعدد الأعمدة يتوافق مع عدد هوائيات الإرسال M_T ، البعد الجانبي يتوافق مع M_R .

يمكن بسهولة تكييف الصيغة (1) لحالات MIMO الخاصة.


MISO (خرج M intiple S ingle Output - عدة هوائيات للإرسال وواحدة استقبال):


y = \ sqrt {\ frac {P} {M_T}} \ mathbf {h} \ mathbf {s} + n \ qquad (4)

حيث \ mathbf {h} هو ناقل 1 \ مرات M_T .


SIMO ( S ingle Single M ultiple O utput - عدة هوائيات استقبال وهوائي إرسال واحد):


y = \ sqrt {P} \ mathbf {h} s + \ mathbf {n} \ qquad (5)

حيث \ mathbf {h} هو ناقل M_R \ المرات 1


SISO ( S ingle Single S ingle O utput - هوائي واحد في جانبي الاستقبال والإرسال):


y = \ sqrt {P} hs + n \ qquad (6)

يبدو أن تكون بسيطة.


يمكن تقسيم كل الاعتبارات الإضافية إلى حالتين كبيرتين: معلومات حالة القناة (CSI - معلومات حالة القناة) غير معروفة للمرسل ( CU - قناة غير معروفة) وتعرف معلومات حالة القناة إلى المرسل ( CK - Channel معروف).


أعلاه ، قمنا بفحص الحالة عندما تكون القناة غير معروفة للمرسل (علبة الحلقة المفتوحة ، الإرسال بدون تعليق). بمعنى آخر ، نظرًا لعدم توفر المعلومات اللازمة ، لا يمكننا اختيار أي اتجاه فعال ، وبالتالي نتبع أبسط مسار: نرسل القدرة المتساوية عبر جميع الهوائيات (المسارات ، مسارات الانتشار). لذلك ، كسب كل كسب مسار هو 1 :


\ gamma_i = 1 ، \ quad i = 1،2 ، .. M_T \ qquad (7)

ومع ذلك ، نكرر: نحن نريد تخصيص المزيد من القدرة للقنوات الجيدة (مسارات الانتشار) وخفض الطاقة للقنوات السيئة .


السؤال الذي يطرح نفسه: كيف يمكن توزيع السلطة بشكل فعال؟


إذا كانت القناة معروفة ( حالة حلقة مغلقة - مع تعليقات) ، فيمكننا استخدام سيناريوهات نقل متقدمة مع بعض خوارزميات معالجة الإشارات الإضافية. على سبيل المثال ، مع النهج الخطية مثل الترميز المسبق والمعالجة اللاحقة .


سوف نفهم معنى المصطلحين الأخيرين.


إذا كان لدينا CSI على جانب الإرسال ، أي قالب \ mathbf {H} ، يمكننا معالجة هذه المصفوفة حسابيا. على سبيل المثال ، تطبيق خوارزمية SVD (تحليل قيمة المفرد).



لاحظ أن المصفوفة \ mathbf {\ Sigma} هي مصفوفة قطرية ، وعناصر قطريها (القيم الفردية) هي ، في جوهرها ، معاملات الإرسال لمسارات الانتشار الفريدة. بمعنى آخر ، إذا حققنا مضاعفة الإشارة لدينا بمصفوفة من القيم الفردية \ mathbf {\ Sigma} بدلا من قناة كاملة \ mathbf {H} ، ثم تتحول قناة MIMO إلى مجموعة من قنوات SISO المتوازية.


لذلك يجب أن تكون مصفوفة الترميز الخطي (مرشح) \ mathbf {F} = \ mathbf {V} _s ، ومصفوفة ما بعد المعالجة الخطية (المستخلص) \ mathbf {D} = \ mathbf {U} ^ H_s ( H لتقف على الاقتران الهرمي).


من الواضح ، بالنسبة للقضية مع قناة غير معروفة \ mathbf {F} و \ mathbf {D} مصفوفات الهوية المتساوية.

الآن ، مع العلم بكل ما سبق ، دعونا نعيد تعريف نموذج الإشارة المستلمة:


\ mathbf {Dy} = \ mathbf {D} \ left (\ sqrt {\ frac {P} {M_t}} \ mathbf {H} \ mathbf {F} \ mathbf {s} + \ mathbf {s} + \ mathbf {n} \ right) = \ mathbf {U} ^ H_s \ mathbf {y} = \ sqrt {\ frac {P} {M_t}} \ mathbf {U} ^ H_s \ mathbf {H} \ mathbf {V} _s \ mathbf {s} + \ mathbf {U} ^ H_s \ mathbf {n} = \ sqrt {\ frac {P} {M_t}} \ mathbf {\ Sigma} _s \ mathbf {s} + \ mathbf {\ hat {n}} = \ mathbf {\ hat {y}} \ qquad (8)

لاحظ أن:


  • \ mathbf {\ hat {n}} لديه نفس الخصائص الإحصائية مثل \ mathbf {n} .
  • قيم إيجن \ mathbf {HH} ^ H هي المربعات للقيم الفردية لمصفوفة القناة \ mathbf {H} ( \ sigma_i = \ sqrt {\ lambda_i} ).

بشكل تخطيطي ، يمكن تمثيل هذا كـ:


شوب

التين. 1. مخطط الترميز المسبق وما بعد المعالجة [1 ، صفحة 67].


parall

التين. 2. مخطط التحلل مشروط \ mathbf {H} عندما تكون القناة معروفة للمرسل والمستقبل [1 ، صفحة 67].


يتم تفكيك الأساسيات - يمكننا المتابعة مباشرةً إلى النطاق الترددي !


الإنتاجية (القدرات)


أعتقد أن كل من درس نظرية المعلومات يتذكر أن مصطلح النطاق الترددي جاء إلينا من هذا النظام المعين. عادةً (في عمر الطالب) ، تم التركيز على الحالة الكلاسيكية لقناة AWGN ، ولكن يمكن اشتقاق الصيغة نسبيًا لحالة قناة MIMO مع الخبو.


لكي لا أعد كتابة الحسابات من الكتب مرة أخرى ، حاولت أن أرتب كل شيء بشكل أو بآخر ، ومن ناحية - لإعطاء صيغ الحياة ، إذا جاز التعبير. آمل أن يكون هذا التنسيق أقل مملة.

لذلك ، مرة أخرى ، نكتب نموذج الإشارة:



الآن دعنا ننتقل إلى تحديد الإنتاجية من خلال المعلومات المتبادلة .



نكتب مصفوفة التوفيق التلقائي للإشارة المستلمة ومكوناتها:



ونحن نستخدمها في تحديد الانتروبيا التفاضلية :



البديل (4) و (5) في (2):



والآن نستبدل (6) في (1):



نواصل المنطق. خذ الحالة الأولى: القناة غير معروفة ( C hannel U nknown). هذا يعني بالنسبة لنا أنه من المستحيل اختيار الاتجاه الأمثل للإرسال ، وبالتالي فإن الإشارات المرسلة ستكون مستقلة وستكون لها نفس القدرة (مدعومة بشكل متساوٍ). بناءً على الحالة القصوى ( Tr \ {\ mathbf {R} _ {ss} \} = M_T ) ، يمكننا أن نأخذ مصفوفة التوفيق التلقائي للأحرف المرسلة مساوية لمصفوفة الهوية. ثم لدينا:



نستخدم خاصية المحددات التالية:


det \ left (\ mathbf {I} _m + \ mathbf {AB} \ right) = det \ left (\ mathbf {I} _n + \ mathbf {BA} \ right)

هذه هي حالتنا ، ويمكننا تبديل المصفوفات بحيث \ mathbf {Q} \ mathbf {Q} ^ H = \ mathbf {I} (من خصائص EVD). ستبقى:



بالانتقال من المصفوفات إلى المبالغ ، لدينا:



توضح هذه الصيغة مرة أخرى النهج المتبع في اعتبار MIMO كقنوات SISO موازية.
بالنسبة إلى القناة المعروفة ( C hannel K nown) ، ستتم إضافة أوزان الهوائي إلى الصيغة:



نكتب أيضًا صيغًا للحالات الخاصة:



ملاحظة :
بالنسبة لحالات SIMO و MISO ، لم تذهب سدى إلى ظهور مربعات قاعدة Frobenius في السجل || \ mathbf {h} || _F ^ 2 - من وجهة نظر رياضية ، فهي تعادل القيم الذاتية \ mathbf {h} \ mathbf {h} ^ H . لذلك ، إذا كنت بحاجة إلى حساب شيء يدويًا بسرعة - فإليك طريقة.

حسنًا ، آمل ألا تتداخل خطي في اللغة الإنجليزية ولغتي الإنجليزية كثيراً مع إدراك المعلومات ، لكن مع ذلك ، دعنا نتحدث عن النقطة الرئيسية :


  • نعم ، يمكن اعتبار عرض نطاق قناة MIMO هو مجموع عرض نطاق قناة SISO .
  • ومع ذلك ، يقتصر هذا المبلغ حسب رتبة القناة!

خوارزمية صب الماء


كما يتضح من صيغة النطاق الترددي المعروفة على جانب الإرسال للقناة (CK - القناة المعروفة) ، يمكن تحسين توزيع الطاقة عبر الهوائيات. للقيام بذلك ، نستخدم خوارزمية صب الماء ( ملء بالماء ) [1 ، ص. 68-69]:


import numpy as np from numpy import linalg as LA import matplotlib.pyplot as plt def waterpouring(Mt, SNR_dB, H_chan): SNR = 10**(SNR_dB/10) r = LA.matrix_rank(H_chan) H_sq = np.dot(H_chan,np.matrix(H_chan, dtype=complex).H) lambdas = LA.eigvals(H_sq) lambdas = np.sort(lambdas)[::-1] p = 1; gammas = np.zeros((r,1)) flag = True while flag == True: lambdas_r_p_1 = lambdas[0:(r-p+1)] inv_lambdas_sum = np.sum(1/lambdas_r_p_1) mu = ( Mt / (r - p + 1) ) * ( 1 + (1/SNR) * inv_lambdas_sum) for idx, item in enumerate(lambdas_r_p_1): gammas[idx] = mu - (Mt/(SNR*item)) if gammas[rp] < 0: #due to Python starts from 0 gammas[rp] = 0 #due to Python starts from 0 p = p + 1 else: flag = False res = [] for gamma in gammas: res.append(float(gamma)) return np.array(res) 

الاختبار:


 Mt = 3 SNR_db = 10 H_chan = np.array([[1,0,2],[0,1,0], [0,1,0]], dtype = float) gammas = waterpouring(Mt, SNR_db, H_chan) print('Rank of the matrix: '+str(LA.matrix_rank(H_chan))) print('Gammas:\n'+str(gammas)) >>> Rank of the matrix: 2 >>> Gammas: >>> [1.545 1.455] 

حسنا ، يبدو معقولا:
1) عدد هوائيات الإرسال المعنية يساوي رتبة القناة ؛
2) مجموع أوزان الهوائيات يساوي عدد هوائيات الإرسال.


حالتان تحد


والآن ، دعونا نتشتت قليلاً ونحل مشاكل التفاهم.


دعنا نجد ، على سبيل المثال ، ما هي المعاملات التي ستكون مساوية لها \ gamma_i مع SNR تميل إلى + \ infty و - \ infty (على نطاق لوغاريتمي ، بالطبع ، لأنه لا توجد قوى سلبية).


نذكر صيغة المراسلات بين الديسيبل والأوقات:


SNR_ {dB} = 10log_ {10} \ left (\ frac {S} {N} \ right)

حيث S - قوة الإشارة المرسلة (بالنسبة لمهامنا ، فهي تعادل طاقة الرمز E_s ) ، و N - قدرة الضوضاء (في مشكلتنا تساوي الكثافة الطيفية للضوضاء N_0 ).


لذلك على مقياس خطي سيكون:


\ frac {E_s} {N_0} \ equiv \ frac {S} {N} = 10 ^ {SNR_ {dB} / 10}

نحن ننظر إلى الصيغ الأساسية للخوارزمية:


\ mu = \ frac {M_T} {(r-p + 1)} \ left [1 + \ frac {N_0} {E_s} \ sum_ {i = 1} ^ {r-p + 1} \ frac {1} {\ lambda_i} \ اليمين]

حيث ص هو التكرار بدءا من 1 ، ص هو رتبة مصفوفة القناة ، \ lambda_i - القيمة الذاتية الأولية لـ "مربع" مصفوفة القناة. يتم حساب Gammas وفقًا للصيغة التالية:


\ gamma_i = \ left (\ mu - \ frac {M_TN_0} {E_s \ lambda_i} \ right) \ quad i = 1،2، ...، r-p + 1

نبدأ في التفكير:


إذا SNR_ {dB} \ إلى + \ infty ، ثم \ frac {E_s} {N_0} \ إلى + \ infty . لذلك، \ frac {N_0} {E_s} \ إلى 0 . بالنسبة للتكرار الأول ، يبقى:


\ mu = \ frac {M_T} {r}

بديلاً عن جاما:


\ gamma_i = \ mu = \ frac {M_T} {r} \ quad i = 1،2 ، ... ، r

نحن نلخص:


مع طاقة نقل كبيرة بلا حدود أو ضوضاء متناهية الصغر ، لا يجب اختراع أي شيء خاص ، دعنا نقول ذلك - نحن نوزع بالتساوي القوة بين هوائيات الإرسال (مع مراعاة مصفوفة القناة).


نحن سبب إضافي:


وماذا عن حالة SNR تميل إلى - \ infty ؟ هنا لن نذهب إلى الرياضيات ، بل سنقوم بالسبب المنطقي: هذه الحالة تقابل إما ضوضاء كبيرة بلا حدود أو قدرة إرسال صفرية. لذلك ، بهذه الطريقة وهذا ، نظامنا ، النظر ، لا يعمل. لذلك ، يختفي السؤال مع gammas تلقائيًا ...


هذه هي الأسئلة التي تظهر في بعض الأحيان في امتحان مع أستاذ.


حساب الإنتاجية (أخيرًا!)


 def siso_capacity(H_chan, SNR_dB): SNR = 10**(SNR_dB/10) c = np.log2(1 + SNR*(np.abs(H_chan)**2)) return c def openloop_capacity(H_chan, SNR_dB): SNR = 10**(SNR_dB/10) Mt = np.shape(H_chan)[1] H_sq = np.dot(H_chan,np.matrix(H_chan, dtype=complex).H) lambdas = LA.eigvals(H_sq) lambdas = np.sort(lambdas)[::-1] c = 0 for eig in lambdas: c = c + np.log2(1 + SNR*eig/Mt) return np.real(c) def closedloop_capacity(H_chan, SNR_dB): SNR = 10**(SNR_dB/10) Mt = np.shape(H_chan)[1] H_sq = np.dot(H_chan,np.matrix(H_chan, dtype=complex).H) lambdas = LA.eigvals(H_sq) lambdas = np.real(np.sort(lambdas))[::-1] c = 0 gammas = waterpouring(Mt, SNR_dB, H_chan) for idx, item in enumerate(lambdas): c = c + np.log2(1 + SNR*item*gammas[idx]/Mt) return np.real(c) Mr = 4 Mt = 4 H_chan = (np.random.randn(Mr,Mt) \ + 1j*np.random.randn(Mr, Mt))/np.sqrt(2) #Rayleigh flat fading c = openloop_capacity(H_chan, 10) print(c) c = closedloop_capacity(H_chan, 10) print(c) c = siso_capacity(H_chan[0,0], 10) print(c) >>> 11.978909197556913 >>> 12.342571770086721 >>> 3.9058582578551193 

يبدو أن العمل. ننتقل إلى المزيد من التقييمات الموضوعية.


قدرة ارجوديك


كما يتضح من الأمثلة أعلاه ، نحن نعمل مع عمليات عشوائية. وبصراحة ، من الخطأ استخلاص أي استنتاجات حول العمليات العشوائية في تطبيق واحد. حتى إذا كانت القناة ثابتة بالمعنى الإحصائي ، فمن الضروري إجراء بعض التوسعات على مجموعة كبيرة بدرجة كافية.


هنا مفهوم القدرة ergodic مفيد لنا:


\ hat {C} = E \ left \ {C \ right \} \ qquad (9)

حيث E \ {* \} يدل على حصيرة. التوقع (القيمة المتوقعة).


نحن النمذجة.


 Mr = 4 Mt = 4 counter = 1000 SNR_dBs = [i for i in range(1, 21)] C_MIMO_CU = np.empty((len(SNR_dBs), counter)) C_MIMO_CK = np.empty((len(SNR_dBs), counter)) C_SISO = np.empty((len(SNR_dBs), counter)) C_SIMO = np.empty((len(SNR_dBs), counter)) C_MISO_CU = np.empty((len(SNR_dBs), counter)) C_MISO_CK = np.empty((len(SNR_dBs), counter)) for c in range(counter): H_MIMO = (np.random.randn(Mr,Mt) + 1j*np.random.randn(Mr, Mt))/np.sqrt(2) H_SISO = H_MIMO[0,0] H_SIMO = H_MIMO[:,0].reshape(Mr,1) H_MISO = H_MIMO[0,:].reshape(1,Mt) for idx, SNR_dB in enumerate(SNR_dBs): C_MIMO_CU[idx, c] = openloop_capacity(H_MIMO, SNR_dB) C_MIMO_CK[idx, c] = closedloop_capacity(H_MIMO, SNR_dB) C_SISO[idx, c] = siso_capacity(H_SISO, SNR_dB) C_SIMO[idx, c] = openloop_capacity(H_SIMO, SNR_dB) C_MISO_CU[idx, c] = openloop_capacity(H_MISO, SNR_dB) C_MISO_CK[idx, c] = closedloop_capacity(H_MISO, SNR_dB) C_MIMO_CU_erg = np.mean(C_MIMO_CU, axis=1) C_MIMO_CK_erg = np.mean(C_MIMO_CK, axis=1) C_SISO_erg = np.mean(C_SISO, axis=1) C_SIMO_erg = np.mean(C_SIMO, axis=1) C_MISO_CU_erg = np.mean(C_MISO_CU, axis=1) C_MISO_CK_erg = np.mean(C_MISO_CK, axis=1) plt.figure(figsize=(7, 5), dpi=600) plt.plot(SNR_dBs, C_MIMO_CU_erg,'g-o', label='$M_R=4$, $M_T=4$ (CU)') plt.plot(SNR_dBs, C_MIMO_CK_erg,'g-v', label='$M_R=4$, $M_T=4$ (CK)') plt.plot(SNR_dBs, C_MISO_CU_erg, 'm-o', label='$M_R=1$, $M_T=4$ (CU)') plt.plot(SNR_dBs, C_MISO_CK_erg, 'm-v', label='$M_R=1$, $M_T=4$ (CK)') plt.plot(SNR_dBs, C_SISO_erg, 'k-', label='$M_R=1$, $M_T=1$') plt.plot(SNR_dBs, C_SIMO_erg, 'c-', label='$M_R=4$, $M_T=1$') plt.title("Ergodic Capacity") plt.xlabel('SNR (dB)') plt.ylabel('Capacity (bps/Hz)') plt.legend() plt.minorticks_on() plt.grid(which='major') plt.grid(which='minor', linestyle=':') plt.show() 


الشكل 3. منحنيات عرض النطاق الترددي لأنظمة نقل مختلفة. قارن مع [1 ، ص. 74].


لذلك نحن نرى ذلك


  • من المتوقع أن تكون حالة MIMO أعلى من غيرها ، ومع زيادة نسبة الإشارة إلى الضوضاء (SNR) ، تقل الحاجة إلى معرفة مصفوفة القناة (انظر المثال مع اللانهاية).
  • SIMO أفضل من MISO إذا كان المرسل لا يعرف القناة (يتم تقاسم القدرة في MISO عبر جميع الهوائيات ، ولكن ليس على النحو الأمثل) ويتزامن مع MISO في حالة قناة معروفة .
  • SISO ومن المتوقع أن جلد في الذيل.

ويحكم جلالة الملك على رتبة مصفوفة القناة ، والذي لا يسمح لنا بمقارنة الزيادة في عدد الهوائيات بشكل لا لبس فيه مع زيادة سرعة الإرسال.


مثل هذه الأشياء.


أدب


(على الرغم من وجود كتاب واحد ، ولكن ما كتاب!)


  1. بولراج ، أروجواسوامي ، روهيت نبار ، ودانانجاي غور.
    مقدمة في الاتصالات اللاسلكية الزمكان. مطبعة جامعة كامبريدج ، 2003.

Source: https://habr.com/ru/post/ar448570/


All Articles