[DELTA67]

السلام عليكم
ملاحظة: هذا الموضوع يخص الحاسوب لكن يمكن استغلال المعلومات لاستعمالها مع المكروكنترولر.

تحويل فورييه السريع Fast Fourier Transform أو اختصارا FFT هو مجموعةُ خوارزميات تُستعمَلُ لحساب تحويل فورييه المقطَّع Discrete Fourier Transform أو ( DFT ) لدالّة (قيم الاشارة المقطعة بعد تحويلها من تماثلية الى رقمية باستعمال تردد معين) و ذلك بسرعة كبيرة ( نتيجة تخفيض العمليات الحسابية اللازمة للمهمة بمئات المرات) و بالتالي يتم التحويل في وقت قصير جدا.
تستخدم هذه الطريقة لتحويل دالّة رياضية بمتغير حقيقي وذات قيم مركّبة في نطاق الزمن إلى دالّة أخرى مكافئة في نطاق التّردّد. و بالتالي فهي تستعمل كثيرا لاجراء التحليل الطيفي لاشارة و هذا ما سنراه ان شاء الله.

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

هذه بعض النقاط المهمة التي يجب أخذها بعين الاعتبار لاستعمال الـ FFT لهذه المهمة:

- نستعمل في كل مرة عينة sample من الاشارة طولها n حيث n =2^N أي أن طول العينة التي نرسلها للتحويل يجب أن يكون مساويا لقوة صحيحة للرقم 2 . مثلا 64 , 128, 256, 512, 1024 .... كلما كان طول العينة أكبر زادت دقة التحويل لكن بالمقابل يزيد الزمن اللازم لاجرائه.
- بما أننا نتعامل مع اشارة (دالة) لمتغير حقيقي فاننا نعوض الاجزاء التخيلية لاشارتنا بأصفار. كيف ذلك؟ توضع العينة في شعاع (جدول) Array القيمة الاولى له هو أول قيمة لعينتنا (الجزء الحقيقي) تكون متبوعة بصفر (الجزءالتخيلي) ثم نضع في القيمة الثالثة للشعاع القيمة الثانية لاشارتنا متبوعة بصفر ... و هكذا. في النهاية نتحصل على شعاع طوله ضعف العينة المراد تحليلها.

- ناتج عملية التحويل (الخرج) سيخزن في نفس الشعاع (الجدول) الذي استعملناه للدخل( أي لامداد تحويل فورييه بقيم عينة اشارتنا). يحتوي هذا الشعاع على قيم الترددات المكونة للاشارة, كل تردد يكون على جزئين حقيقي و تخيلي.
- أول قيمة للخرج (أي التردد عند الدليل 0 = i في الشعاع) هي قيمة خاصة حيث تمثل المتوسط لكل قيم ترددات العينة. القيم الأخرى من الدليل i=1 الى الدليل i=n/2 تمثل قيم الترددات المشكلة للاشارة و التي تحسب بالعلاقة:

f = samplingRate * i / n تكون قيمة التردد بالهرتز.

لاحظ أن الترددات التي يمكننا معرفتها باستعمال تحويل فورييه تعتمد على طول العينة n و قيمة التردد المستعمل لترقيم الاشارة samplingRate. فمثلا اذا كانت قيمة هذا الأخير 44100 Hz و طول عينتنا 1024 =n فان التردد عند الدليل i=1 هو 44100 Hz * 1 / 1024 = 43.07 Hz بقية الترددات هي مضاعفات لهذا التردد فمثلا قيمة التردد عند الدليل i = 17 هي 43.07*17 = 732.13 Hz.
التردد عند الدليل i=n/2 هو حالة خاصة أيضا حيث يمثل تردد Nyquist.
- لحساب طويلة magnitude و عمدة angle التردد عند الدليل i (من شعاع الخرج [...]A) نستعمل العلاقات التالية:

      double magnitude = sqrt (A[i]*A[i]+A[i+1]*A[i+1]);

    double angle = atan2 (A[i], A[i+1] );

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

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

1- توليد الاشارة: يمكن استعمال اشارة جيبية أو مربعة ذات تردد محدد و مرقمة
عند تواتر معين (44100 هرتز) و هو المستعمل غالبا للانتاج الصوتي.

اشارة جيبية: ذات تردد واحد .

int i;	
float sample, amplitude, A[2048];

  /* géneration d'un signal sinusoidal 1000hz echantilloné à 44100*/


	for (i=0; i<1024; i++)
	{ sample= 32000*sin((i/44.1)*2*PI); // 44100:1000=44.1 (44.1 echantillons par periode)
	  A[2*i]=sample;
	  A[2*i+1]=0;
	}

لاحظ أن طول العينة 1024 و بالتالي طول الشعاع الذي نخزنها به هو 2048 , كل قيمة من العينة تحتل خانتين من الشعاع (خانة للجزء الحقيقي أي قيمة العينة ذاتها و الخانة التي تليها مباشرة للجزء التخيلي و الذي قلنا سابقا أننا نضعه صفرا)

اشارة مربعة: تحتوي على تردد أصلي Fundamental و عدة توافقيات Harmonics .

/* *******  ou un signal  rectangle  (plusieurs harmoniques) */ 

int i, j=5;
float sample, amplitude, A[2048];

  for( i = 0 ; i < 1024 ; i++)
    {
      if( i % 44 == 0 ) // 44 echantillons par période
	{
	  if ( j == 5 ) // amplitude = 5 volt
	    j = -5;
	  else
	    j = 5;
	}
      A[2*i] = j;
	  A[2*i+1] =0;
   }

ملاحظة: يمكن استخدام قطعتي الكود لتوليد اشارة جيبية أو مربعة لكتابة برنامج مولد اشارات ذات تردد منخفض . يمكن اضافة الهيدر Header الخاص لتكوين ملف نوع WAV .

2- نقوم بمعالجة اشارتنا باستعمال الـ FFT:
نقوم باستدعاء الغوريتم التحويل بواسطة الأمر:

FFT(A,1024,1);

حيث FFT هو اسم الالغوريتم و A هو الشعاع الذي يحتوي على عينة اشارتنا و 1024 طولها و الواحد 1 يعلم الالغوريتم على اجراء تحويل فورييه ( نضع ناقص واحد اذا كنا نريد التحويل العكسي).
هذا هو الالغوريتم المستعمل و هو مأخوذ من كتاب Numerical Recipes In C مع بعض التعديل و التحسين من طرف أحد الأشخاص على النت.

void FFT(float data[], unsigned long N, int isign)
{

	//data -> float array that represent the array of complex samples
//N -> number of samples (N^2 order number) 
//isign -> 1 to calculate FFT and -1 to calculate Reverse FFT
	
	
	//variables for trigonometric recurrences
    unsigned long n,mmax,m,j,istep,i;
    double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
//the complex array is real+complex so the array 
    //as a size n = 2* number of complex samples
    // real part is the data[index] and 
    //the complex part is the data[index+1]
    n=N * 2; 

    //binary inversion (note that the indexes 
    //start from 0 witch means that the
    //real part of the complex is on the even-indexes 
    //and the complex part is on the odd-indexes
    j=0;
    for (i=0;i<n/2;i+=2) {
        if (j > i) {
            //swap the real part
            SWAP(data[j],data[i]);
            //swap the complex part
            SWAP(data[j+1],data[i+1]);
            // checks if the changes occurs in the first half
            // and use the mirrored effect on the second half
            if((j/2)<(n/4)){
                //swap the real part
                SWAP(data[(n-(i+2))],data[(n-(j+2))]);
                //swap the complex part
                SWAP(data[(n-(i+2))+1],data[(n-(j+2))+1]);
            }
        }
        m=n/2;
        while (m >= 2 && j >= m) {
            j -= m;
            m = m/2;
        }
        j += m;
    }

//Danielson-Lanzcos routine 
    mmax=2;
    //external loop
    while (n > mmax)
    {
        istep = mmax<<  1;
        theta=isign*(6.28318530717959/mmax);
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        //internal loops
        for (m=1;m<mmax;m+=2) {
            for (i= m;i<=n;i+=istep) {
                j=i+mmax;
                tempr=wr*data[j-1]-wi*data[j];
                tempi=wr*data[j]+wi*data[j-1];
                data[j-1]=data[i-1]-tempr;
                data[j]=data[i]-tempi;
                data[i-1] += tempr;
                data[i] += tempi;
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }
}

- الآن أصبح الشعاع [..] A يحتوي على الترددات المكونة لاشارتنا ما بقي لنا سوى حساب طويلة كل منها و عرضه على الشاشة . لا يتم عرض قيمة التردد بل دليله و يمكن حساب التردد كما رأينا سابقا:

/* Affiche la TF du signal d'entrée */

  i = 0 ; 
  while ( i < 256) //les 256 premiers amplitudes
    {   
        amplitude = sqrt(A[i]*A[i]+A[i+1]*A[i+1]);
		printf(" l'amplitude à l'indice %ld est: %f \n",i, amplitude);
		i=i+2;
    }

من خلال القيم المعروضة نعرف الترددات التي تحتوي عليها اشارتنا و طويلة كل منها.

-----------------
المراجع:

1- أعتمدت على هذه الصفحة:

http://www.codeproject.com/Articles/...-FFT-algorithm

2- كتاب Numerical Recipes In C يمكنكم تصفحه مجانا على النت هنا:

http://www.nr.com

كما يمكن بالبحث الدؤوب ايجاد نسخ منه يمكن تحميلها.

أحبــــــــــك و الله يا مرســـــــــــــي

--------------------------------------------------------------------------------------
مـــــواضـــيعــي: (أنقـــر على العنوان لتصفح الموضوع)

-هيا نصنع دارة محول صوت مونو الى "شبه ستيريو"Mono to Pseudo STEREO.

- كيف تستعمل شاشة تلفون نوكيا 3310 أو 3410؟

- ما رأيكم في مشروع يناء حاسوب موافق للـ IBM PC 5150 ؟؟

- مبرمجة الـــ PIC داخل الدارة In Circuit.

- وصل بطاقة SD أو MMC بالـــ ATMEGA8.

- مبرمجة لعائلة الــ AVR سهلة جدا جدا جدا !!!

- دارة بسيطة جدا لعرض نص على شاشة التلفزيون.

- مبرمجات PIC بسيطة جدا.

- أشعل شمعة الكترونية، استرخي تحت وقع زخات المطرثم نم نوما هادئا!!!

- مبرمجـــة الــ ATMEL AVRs عن طريق الــ USB .

- أضف شاشة عرض LCD الى مشاريعك.

- ما رأيكم في مشروع دايزك DISEQC ؟؟

- هل تعرف LTSPICE IV ؟ برنامج محاكاة احترافي قوي و مجاني.

- تعالوا نتعرف على الصوت المحيطي SURROUND SOUND .

- مجموعة من الحيل "العفسات" للمحترفين و الهواة.

- مدخل لاستخدام البورت USB. أستعد لنهاية البورتات LPT و RS232.

- كيف تصنع جهاز استقبال راديو بترانزستور واحد؟؟

- كيفية استعمال تحويل فورييه السريع FFT لايجاد طيف اشارة.

[DELTA67]

و هذا هو الكود كاملا:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.1415926535897932385
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void FFT(float data[],unsigned long N,int isign);


int main()
{ 
	int i, j=5; // j pour le signal carre
	float sample, amplitude, A[2048];



  /* géneration d'un signal sinusoidal 1000hz echantilloné à 44100*/


	for (i=0; i<1024; i++)
	{ sample= 32000*sin((i/44.1)*2*PI); // 44100:1000=44.1 (44.1 echantillons par periode)
	  A[2*i]=sample;
	  A[2*i+1]=0;
	}
 

	 /* *******  ou un signal  rectangle  (plusieurs harmoniques) 



  for( i = 0 ; i < 1024 ; i++)
    {
      if( i % 44 == 0 ) // 44 echantillons par période
	{
	  if ( j == 5 ) // amplitude = 5 volt
	    j = -5;
	  else
	    j = 5;
	}
      A[2*i] = j;
	  A[2*i+1] =0;

    }     *******  */



FFT(A,1024,1);  


/* Affiche la TF du signal d'entrée */

  i = 0 ; 
  while ( i < 256) //les 256 premiers amplitudes
    {   
        amplitude = sqrt(A[i]*A[i]+A[i+1]*A[i+1]);
		printf(" l'amplitude à l'indice %ld est: %f \n",i, amplitude);
		i=i+2;
    }
 
system("PAUSE");
return 0;
}
 


void FFT(float data[], unsigned long N, int isign)
{

	//data -> float array that represent the array of complex samples
//N -> number of samples (N^2 order number) 
//isign -> 1 to calculate FFT and -1 to calculate Reverse FFT
	
	
	//variables for trigonometric recurrences
    unsigned long n,mmax,m,j,istep,i;
    double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
//the complex array is real+complex so the array 
    //as a size n = 2* number of complex samples
    // real part is the data[index] and 
    //the complex part is the data[index+1]
    n=N * 2; 

    //binary inversion (note that the indexes 
    //start from 0 witch means that the
    //real part of the complex is on the even-indexes 
    //and the complex part is on the odd-indexes
    j=0;
    for (i=0;i<n/2;i+=2) {
        if (j > i) {
            //swap the real part
            SWAP(data[j],data[i]);
            //swap the complex part
            SWAP(data[j+1],data[i+1]);
            // checks if the changes occurs in the first half
            // and use the mirrored effect on the second half
            if((j/2)<(n/4)){
                //swap the real part
                SWAP(data[(n-(i+2))],data[(n-(j+2))]);
                //swap the complex part
                SWAP(data[(n-(i+2))+1],data[(n-(j+2))+1]);
            }
        }
        m=n/2;
        while (m >= 2 && j >= m) {
            j -= m;
            m = m/2;
        }
        j += m;
    }

//Danielson-Lanzcos routine 
    mmax=2;
    //external loop
    while (n > mmax)
    {
        istep = mmax<<  1;
        theta=isign*(6.28318530717959/mmax);
        wtemp=sin(0.5*theta);
        wpr = -2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        //internal loops
        for (m=1;m<mmax;m+=2) {
            for (i= m;i<=n;i+=istep) {
                j=i+mmax;
                tempr=wr*data[j-1]-wi*data[j];
                tempi=wr*data[j]+wi*data[j-1];
                data[j-1]=data[i-1]-tempr;
                data[j]=data[i]-tempi;
                data[i-1] += tempr;
                data[i] += tempi;
            }
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }
}

[سعيد قادر]

جزاك الله كل الخير يا بشمهندس هذا الموضوع من المواضيع المحببة لدى بارك الله فيك و شكرا جزيلا لك

[emad_e.m.s]

بارك الله فيك الله يعطيك العافية
هل من الضروري أن يحتوي المايكروكونترولر على وحدة معالجة رقمية(DSP) إذا أردنا تطبيق هذه العملية؟؟

[DELTA67]

اقتباس المشاركة الأصلية كتبت بواسطة emad_e.m.s بارك الله فيك الله يعطيك العافية هل من الضروري أن يحتوي المايكروكونترولر على وحدة معالجة رقمية(DSP) إذا أردنا تطبيق هذه العملية؟؟

و فيكم بارك أخوتي.

لا ليس من الضروري أن يحتوي المايكروكونترولر على وحدة معالجة رقمية(DSP) لكن استطاعته القيام بعمليات الضرب مباشرة هذا سيسرع من الحسابات طبعا.

مثلا يمكن اجراء الـ FFT بمكروكنترولر نوع AVR عادي كما تراه هنا:

http://elm-chan.org/works/akilcd/report_e.html

أو هنا:

[emad_e.m.s]

مشكور على هذه الجهود الله يعطيك العافية

[ennng]

بارك الله فيك اخي الكريم على الطرح المميز عندي بعض الاضافات وطلب اذا سمحت لي.
1- بالنسبة sampling frequency لابد ان يكون ضعف تردد الاشارة الاصلية او المراد معالجتها على سبيل المثال عندنا اشارة ب 10KHz فلابد ان يكون sample ب 20KHz او اكثر ليتحقق تردد Nyquist. اذا كان اقل راح يعمل تشوة في الاشارة الخارجة.
2- FFT مهم جدا في معرفة اشارة noise على سبيل المثال عندك اشارة (صوت او غير ذالك..) وبها ضوضاء انت لاتعرف تردد الضوضاء بالضبط فتعمل FFT لتحديد تردد الاشارة الاصلية من اشارة الضوضاء وبالتالي تستطيع عمل فلتر لكي تتخلص من الضوضاء.
اما طلبي اذا عندك علم ب FIR filter انا اعرف الاشياء النظرية عنه لكن اريد كود وتطبيق عنه في pic16 على سبيل المثال low pass filter وجزاك الله خير

[zamalkawi]

لم أقرأ الكود بتمعن
ولكن لي سؤالين، هل هذا الكود يعمل في الزمن الحقيقي، أم أن الإشارة يجب أن تكون مسجلة؟
ولو يعمل في الزمن الحقيقي، فهل يعمل بتحويل فوريير في زمن قصير short time fourier transform STFT أم يعمل بطريقة أخرى؟