DFT (Discrete Fourier Transform) ۱ تبدیل Z یا تبدیل فوریه x[n] که به صورت X(z) و ) jω X(e نمایش داده می شوند از لحاظ محاسباتی دو مشکل دارند: ۱ -محاسبه ا نها برای دنباله های با طول بی نهایت ۲ -این تبدیل ها تابع متغیرهای پیوسته یا در واقع فرکانس های غیر قابل شمارش هستند. MATLAB متغیرهای پیوسته هستند و غیر قابل شمارش هستند. برای استفاده از z و ω باید دنباله ها را قطع کنیم تا طول ا نها محدود شود و سپس این تبدیل ها را در فرکانس های محدودی محاسبه کنیم. تبدیل DTFT و تبدیل Z از لحاظ عددی قابل محاسبه نیستند. تبدیل DFT نمونه برداری از DTFT در حوزه فرکانس است. X(e jω ) = (x[n]) (x[n]) = X[k] = (x[n]) ω= ۲πk در ابتدا از دنباله های متناوب شروع می کنیم و سری فوریه ا نها ۲ سری فور یه Series=DFS) (Discrete-Fourier x[n] = x[n + k] سیگنال متناوب با دوره : n, k در این صورت: x[n] = ۱ ۱ k=۰ X[k]e j۲πkn X[k]} { ضرایب سری فوریه نامیده می شوند و داریم: X[k] = ۱ n=۰ x[n]e j۲πkn توجه کنید که خود X[k] هم با دوره متناوب است: X[k] = X[k + ] می گیریم X[k] { x[n]} = ۱ n=۰ :DFS معادله ا نالیز یا W = e j۲π x[n]w nk ۱
x[n] { X[k]} = ۱ ۱ n=۰ X[k]W nk معادله سنتز یا :IDFS *مثال ۵-۱ کتاب صفحه ۱۴۳ را ببینید. معادلات ا نالیز و سنتز را می توان به شکل ماتریسی نوشت.: X = W x x = ۱ W X W ماتریس با فرمول زیر است: که ۱ ۱... ۱ W [W kn ۱ W ۱... W ۱, ۰ k, n ۱] =. ۱ W ( ۱)... W ( ۲)۲ function [Xk]=dfs(xn,) n=0:-1; k=0:-1; W=exp(-j*2*pi/); nk=n'*k; Wnk=W.^nk; Xk=xn*Wnk; W یک ماتریس مربعی است و ماتریس DFS نامیده می شود. تابع زیر معادلات ماتریس بالا را پیاده می کند: مثال: >>x[n]=[0,1,2,3]; >>=4; Xk=dfs(xn,); * تابع idfs را در صفحه ۱۴۵ ببینید. مثال ۵-۲ صفحه ۱۴۵ پالس مستطیلی متناوب تبدیل فوریه این سیگنال sinc متناوب است. مثال: پالس متناوب مستطیلی = ۵ L نقطه ای که از هر ۲۰ نقطه تکرار می شود. >>L=5;=20;k=-/2:/2; >>xn=[ones(1,l),zeros(1,-l)]; >>Xk=dfs(xn,); >>magxk=abs([xk(/2+1:) Xk(1:/2+1)]); >>Subplot(2,2,1);stem(k,magXk);axis([-/2,/2,-0.5,5.5]) >>title('dfs of SQ wave:l=5,=20') ۲
*توجه: تابع X[k] dfs را برای ۱ k ۰ حساب می کند که X[k] هم ۲ رسم کنیم که اولا متناوب است. در مثال بالا می خواهیم X[k] را برای k ۲ X[k] دارد. توجه داریم که dfs نقطه است و ثانیا یک شیفت فرکانس نسبت به + ۱ [ ] [ ] پس: [ ] متناوب است با دوره X = ۲ X ۲ + = X ۲ ۲π(/۲). اما در اندیس های k به = π است زیرا: k = ۲ فرکانس ω = π معادل k = ۲ درMATLAB برابر ۱ + /۲ است. جای ۰ تا ۱ از ۱ تا هستند پس همچنین معادل فرکانس صفر است: X[۰] = X[۰ + ] = X[] ۲ k ۱ به خاطر متناوب بودن هم ارزند با ۲ + k ۱ + ۲ k ۱ در اندیس MATLAB ۲ + ۱ k ۰ k ۲ در اندیس MATLAB ۱ k ۲ + ۱ *شکلهای ۵-۲ در صفحه ۱۴۷ را برای مقادیر مختلف L و ببینید. X[k]* ها نمونه های ) jω X(e (تبدیل فوریه یک تناوب) در حوزه فرکانس در هارمونیک های ۲π است. تبدیل DFT از یک سیگنال غیر متناوب: طبق قضایای DFT می دانیم که اگر x[n] غیر X(e jω ) = متناوب باشد و ) jω X(e تبدیل فوریه ا ن باشد: + k= x[k]e jωk توجه کنید که حدود زیگما نامحدود است. اگر از ) jω X(e روی دایره واحد در ω k = ۲πk نمونه بگیریم: فرکانس های X[k] X(e jω ) ω= ۲πk = X(e j ۲πk ) ۳
x[n] = ۱ ۱ k=۰ X[k]e j ۲πk n = + k= x[n k] در این صورت x[n] را IDFTی X[k] و بالعکس X[k] را DFTی نقطه ای x[n] می نامیم. رابطه بالا می گوید IDFT ی X[k] تکرار x[n] از هر نقطه است. MATLAB فقط x[n] را در بازه ۱ n ۰ به ما می دهد. یعنی ( + ) IDFT{X[k]} = x[n]r [n] = x[n k] R [n] k= [n] R یک پنجره مستطیلی است که در بازه ۱ n ۰ برابر یک و بقیه جاها صفر است. مثال ۵-۵ صفحه :۱۵۱ u[n] x[n] = (۰ ۷) n است از X(z) روی دایره واحد نمونه های متساوی الفاصله با =,۵,۱۰,۲۰ ۵۰ بگیرید و تاثیر ا ن را در حوزه زمان بررسی کنید. حل: ۱ X(z) = ۱ ۰ ۷z ۱ = z z ۰ ۷ z ۰ ۷ X[k] = X(z) z=e j ۲πk k Z *فرمول ساخت تبدیل Z سیگنال x[n] اگر x[n] در بازه ۱ n ۰ محدود باشد از نمونه های DFT ی نقطه ای به صورت زیر است: X(z) = ۱ z ۱ X[k] = ۱ W k z ۱ k=۰ با قرار دادن z = e jω می توان از رابطه بالا ثابت کرد (تمرین): ۱ ( X(e jω ) = X[k]Φ ω ۲πk ) k=۰ ۱ Φ(ω) = sin(ω ۲ ) و = ۱.Φ(۰) توجه کنید این فرمول در واقع درون یابی jω( که ) ۲ sin ω e ۲ تبدیل فوریه از نمونه های تبدیل فوریه است و شبیه فرمول درون یابی سیگنال ا نالوگ از نمونه هایش در حوزه زمان است: x a (t) = + n= x[n](f s (t nt s )) *تفاوتهای این دو فرمول را توضیح دهید. توجه کنید Φ(ω) یک متناوب است. *توجه کنید که مهم است بدانید که محاسبه اندیسها در DFT ی نقطه ای به پیمانه است: x[[n]] x(n mod ) ۴
یا *رابطه بالا در صورتی صحیح است که تکرارهای x[n] با هم تداخل نکنند. *تابع rem(n,) در MATLAB باقیمانده تقسیم n بر را که همان ( n) mod است محاسبه می کند. تابع mod(n,) هم موجود است. تفاوت این دو تابع چیست { X[k] ۰ k ۱ X[k] = DFT{x[n]} = [k] X[k]R = در غیر این صورت ۰ X[k] = ۱ n=۰ x[n]w kn, ۰ k ۱ *توجه کنید X[k] = X[[k]] متناوب است و در تمام k های صحیح از تا + تعریف شده حال ا نکه X[k] فقط برای ۱ k ۰ تعریف شده است. *توابع dfs و idfs را می توان dft و idft نامید. مثال ۵-۶: (تمرین) صفحه ۱۵۷: x[n] پالس مستطیلی ۴ نقطه ای است. تبدیل فوریه ) jω X(e و DFT ی ۴ نقطه ای X[k] را رسم کنید و بگویید چرا X[k] فقط در = ۰ k مقدار دارد مثال ۵-۷: اگر بخواهیم DFT ی ۸ نقطه ای از x[n] بگیریم یعنی ۸ نمونه در فرکانس بگیر یم X[k] = X(e jω ) ω= ۲πk ۸ X[k] = X(z) z=e j ۲πk ۸ = ۸ ۰ k ۸ همانطور که از تي وری DFT می دانیم معادل این است که x[n] از هر ۸ نقطه تکرار شود. پس اگر یک تناوب ا ن را در نظر بگیریم: x[n] = {۱, ۱, ۱, ۱, ۰, ۰, ۰, ۰} یعنی باید ۴ تا صفر به انتهای x[n] اضافه کرد و سپس از ا ن DFTی ۸ نقطه ای گرفت. به این کار الصاق صفر( padding (zero می گویند. *نتیجه بسیار مهم: درون یابی ظریف تر در حوزه ی فرکانس و گرفتن نمونه های بیشتر معادل الصاق صفر در حوزه ی زمان است کد MATLAB (برای = ۸ ): ۵
که در زیر مشاهده می شود. >>X=[1,1,1,1, zeros(1,4)]; =8; x=dft(x,); شکل :۱ DFT برای = ۸ برای اینکه نشان دهیم DFT چند نقطه داشته تعداد نقاط را به صورت اندیس می گذاریم پس [k] X ۸ یعنی ۸ DFT نقطه ای و[ k ] Xیعنی ۱۶ DFT ی ۱۶ نقطه ای. دقت کنید که الصاق صفر ) jω X(e را ظریف تر رسم می کند اما به دقت فرکانسی کمک نمی کند. دلیل این امر ا ن است که الصاق صفر اطلاعات به سیگنال x[n] اضافه نمی کند که ما اطلاعات دقیق تری از فرکانس به دست ا وریم. برای به دست ا وردن اطلاعات بیشتر باید تعداد بیشتری از نقاط داده x[n] را از ا زمایش و مشاهدات عملی به دست ا ورد. دنباله ۵۲πn) x[n] = cos(۰ ۴۸πn) + cos(۰ را در نظر بگیرید: می خواهیم طیف این سیگنال با تعداد محدودی از نمونه ها تعیین کنیم. الف) DFT را برای ۱۰ n ۰ و x[n] رسم کنید. ۶
شکل :۲ برای = ۱۶ n = 0:99; x = cos(0.48*pi*n)+cos(0.52*pi*n); n1 = 0:9; y1 = x(1:10); subplot(2,1,1) stem(n1,y1) title('signal x(n), 0 <= n <= 9') xlabel('n') %10-point DFT Y1 = dfs(y1,10); magy1 = abs(y1(1:1:6)); k1 = 0:5; w1 = 2*pi/10*k1; subplot(2,1,2) stem(w1/pi,magy1) title('samples of DTFT Magnitude') xlabel('frequency in \pi units') %100-point DFT of 10-point signal ب) DFT را برای ۱۰۰ n ۰ و x[n] رسم کنید. ۷
شکل :۳ برای = ۱۲۸ n2 = 0:99; y2 = [x(1:1:10) zeros(1,90)]; subplot(2,1,1) stem(n2,y2) title('signal x(n), 0 <= n <= 9 + 90 zeros') xlabel('n') Y2 =dfs(y2,100); magy2 = abs(y2(1:1:51)); k2 = 0:50; w2 = 2*pi/100*k2; subplot(2,1,2) plot(w2/pi,magy2) title('dtft Magnitude') xlabel('frequency in \pi units') axis([0,1,0,10])%100-point DFT of 100-point signal ۸
subplot(2,1,1); stem(n,x); title('signal x(n), 0 <= n <= 99') xlabel('n') X = dfs(x,100); magx = abs(x(1:1:51)); k = 0:1:50; w = 2*pi/100*k; subplot(2,1,2) plot(w/pi,magx) title('dtft Magnitude') xlabel('frequency in \pi units') axis([0,1,0,60]) الف) از اجرای این کد برای ۱۰ DFT نقطه ای با توجه به شکل ۴ خیلی چیزی دستگیرمان نمی شود. شکل ۴: سیگنال و نمونه های ا ن برای = ۱۰ حال به ۱۰ نقطه ۹۰ تا صفر الصاق می کنیم و این دفعه طیف را با دستور plot رسم می کنیم. شکل ۵ نشان می دهد که این سیگنال یک فرکانس غالب در ω = ۰ ۵π دارد. این جواب صحیح نیست زیرا با توجه به سیگنال اصلی در فرکانس ω ۱ = ۰ ۴۸π ۹
شکل ۵: سیگنال و نمونه های ا ن برای = ۱۰۰ و ω ۲ = ۰ ۵۲π داریم تعداد ۱۰ نقطه x[n] اطلاعات کافی را برای تفکیک فرکانس ها ندارد گرچه ) jω X(e ظریف رسم شده است. این ظرافت بیشتر در رسم غلط بودن اطلاعات را تصیح نمی کند. ب) برای گرفتن اطلاعات بیشتر فرکانسی و دقت فرکانسی بیشتر ۱۰۰ نمونه از x[n] را می گیریم. مطابق کد زیر: >> subplot(2,1,1); stem(n,x); >> title( signal x(n), 0 <= n <= 99'); xlabel('n') >> X = dft(x,100); magx = abs(x(1:1:51)); >> k = 0:1:50; w = 2*pi/100*k; >> subplot(2,1,2); plot(w/pi,magx); title('dtft Magnitude'); ۱۰
>> xlabel('frequency in \pi units') دو قله فرکانسی ۰ ۴۸π و ۰ ۵۲π مشاهده می شود که اطلاعات فرکانسی درست می دهد. خواص :DFT خواص DFT در درس DSP به طور مفصل بحث شده است: این بحث ها را اینجا تکرار نمی کنیم. توجه داشته باشید که اندیس ها در DFT به پیمانه ) تعداد نقاط (DFT حساب می شوند. به طور مثال: DFT{x[[ n]] } = X[[ k]] اگر به طور مثال = ۱۶ باشد داریم : x[[ ۲]] ۱۶ = x[[۱۶ ۲]] ۱۶ = x[۱۴] تمام اندیس ها در بازه ۱ n ۰ قرار می گیرند. یا خاصیت مزدوج شدن به صورت زیر است : DFT{x [n]} = X [[ k]] برای دنباله های حقیقی x[n] X[k] دارای تقارن مزدوج متناوب است یعنی : X[[ k]] = X [k] ۰ k ۱ در این صورت کافی است فقط X[k] برای مقادیر زیر حساب شود: ۰ k ۲ در حالتی که زوج است برای ۱۱
۰ k ۱ ۲ در حالتی که فرد است برای مثال: = ۴ X[[ ۱]] ۴ = X [۱] X[۴ ۱] = X [۱] X[۳] = X [۱] X[[ ۲]] ۴ = X [۲] X[۴ ۲] = X [۲] X [۲] = X[۲] x[n] حقیقی است. مثال: = ۵ X[[ ۱]] ۵ = X [۱] X[۵ ۱] = X [۱] X[۴] = X [۱] X[[ ۲]] ۵ = X [۲] X[۵ ۲] = X [۲] X[۳] = X [۲] در مثال = ۴ [۳]X لازم نیست محاسبه شود و از [۱]X به دست می ا ید. در مثال = ۵ [۳]X و [۴]X لازم نیست محاسبه شوند و از [۱]X و[ X[۲ به دست می ا یند. X[ ۲ حقیقی است. [۰]X همیشه حقیقی است و در حالت زوج ] K = ۲ معادل فرکانس ω = π است. محاسبه کانولوشن از روش : overlap-save توضیح و اثبات های ا ن را در DSP بخوانید. این روش برای محاسبه کانولوشن یک سیگنال با طول زیاد با یک فیلتر FIR استفاده می شود. سیگنال به بلوک های نقطه ای تقسیم می شود ولی اگر طول پاسخ ضربه فیلتر M باشد که M < هر بلوک با بلوک قبل هم پوشانی دارد. به طور مشخص اگر از DFT ی نقطه ای استفاده کنیم خروجی هر بلوک نقطه دارد که ۱ M نقطه اول ا ن غلط است و + ۱ M نقطه ا خر ا ن صحیح است. پس برای انتخاب هر بلوک ورودی باید ۱ M نقطه ا خر بلوک قبل ۱ M نقطه اول بلوک فعلی باشد و بلوک ها باهم ۱ M نقطه هم پوشانی داشته باشند. در این صورت خراب شدن ۱ M نقطه خروجی و دور ریختن ا نها اشکالی ایجاد نمی کند. مثال: ۹ n x[n] = n + ۱, ۰ و ۱} ۰, {۱, = h[n] از روش overlap-save با = ۶ برای محاسبه کانولوشن h[n] y[n] = x[n] استفاده کنید. حل: چون = ۳ M باید هر بلوک با بلوک های قبلی خود دو نقطه (۲ = ۱ M) هم پوشانی داشته باشند. x[n] ۱۰ نقطه ای است. برای شروع بلوک اول باید ۲ تا صفر در ابتدای بلوک اضافه کنیم چون دو نقطه اول خراب می شوند. سه بخش نیاز داریم. ۱۲
x ۱ [n] = {۰, ۰, ۱, ۲, ۳, ۴} x ۲ [n] = {۳, ۴, ۵, ۶, ۷, ۸} x ۳ [n] = {۷, ۸, ۹, ۱۰, ۰, ۰} خروجی های [n] y ۲ [n] y ۱ و [n] y ۳ را با روش DFT حساب می کنیم که از هرکدام دو نقطه اول خراب هستند و ۴ باقی سالم هستند. y ۱ = x ۱ ۶ h[n] = { ۳, ۴, ۱, ۲, ۲, ۲ } y ۲ = x ۲ ۶ h[n] = { ۴, ۴, ۲, ۲, ۲, ۲ } y ۳ = x ۳ ۶ h[n] = {۷, ۸, ۲, ۲, ۹, ۱۰ } y = {۱, ۲, ۲, ۲, ۲, ۲, ۲, ۲, ۲, ۲, ۹, ۱۰} y خروجی کل است. اگر از دستور conv که کانولوشن عادی را محاسبه می کند هم استفاده کنید به همان نتیجه می رسیم. کد تابع ovrlpsav در صفحه ۱۸۶ کانولوشن را به روش overlap-save برای دنباله های طولانی ورودی پیاده می کند. می توان قسمت کانولوشن این کد را اصلاح کرد تا با FFT انجام شود. روش overlap-add برای محاسبه کانولوشن: در این روش x[n] به بلوک های غیر هم پوشان تقسیم می شود اما خروجی ها طبق جمع ا ثار به همان صورت هم پوشان جمع می شوند. توضیح کامل تر را در درس DSP ببینید. محاسبه سریع :DFT الگوریتم های جالبی در این زمینه وجود دارد. کتاب DSP را برای حالت هایی که تعداد نقاط DFT را می توان به حاصل ضرب دو عدد تجزیه کرد ببینید. در این حالت ها داده های ورودی را در یک ماتریس می چینید و یک بار DFT سطری تمام سطرها و یک بار DFT ستونی تمام ستونها حساب می شود. اگر به حاصل ضرب بیش از دو عدد تجزیه شود این الگوریتم به صورت بازگشتی تکرار می شود و تعداد محاسبات از ۲ به ۲ log ۲ کاهش می یابد. کد زیر تخمین زمان محاسبه DFT است. >> max = 2048; fft_time=zeros(1,max); >> for n=1:1:max >> x=rand(1,n); >> t=clock;fft(x);fft_time(n)=etime(clock,t); >> end >> n=[1:1:max]; plot(n,fft_time,'.') >> xlabel('');ylabel('time in Sec.'); title('fft execution times') ۱۳
در این کد بردار های تصادفی با طول های = ۱ و = ۲۰۴۸ تولید می شود و FFT ا نها حساب می شود و زمان محاسبه FFT با طول های مختلف اندازه گیری و رسم می شود. اگر [n] x ۱ طول و ۱ [n] xطول ۲ ۲ داشته باشد و بخواهیم کانولوشن [n] x ۲ و [n] x ۱ را حساب کنیم که طول ا ن ۱ ۲ ۱ + است برای سرعت بیشتر این کانولوشن را می توان با FFT ی نقطه ای پیاده کرد که اولین توان ۲ است که از ۱ ۲ ۱ + بزرگتر است یا با ا ن مساوی است. مثال : ۱ = ۳, ۲ = ۱۲ ۲ + ۱ ۱ = ۱۴ = ۲ ۴ = ۱۶ x ۱ [n] x ۲ [n] = IFFT{FFT{x ۱ [n]} FFT{X ۲ [n]}} تابع fftfilt در MATLAB روش سریع پیاده سازی overlap-add است. ۱۴