Implementation of FFT and IFFT in Matlab
in [Matlab]
|
Prev: ??? Unexpected or internal error encountered in "slOneSimStep". Please report this to The MathWorks if you can cause it to recur.
Next: ODE and external function as input
From: Wayne King on 18 Mar 2010 07:37 "Saurabh " <saurabhs(a)student.ethz.ch> wrote in message <hnt295$al5$1(a)fred.mathworks.com>... > well sorry for missing dt which is 1E-6..... I am new to matlab > > but the thing is that if I just use the code > > Pdft = fft(p_i); > plot (abs(Pdft)) > or > Pdft = fft(p_i); > stem(abs(Pdft)) > > I do not get what is expected so I have not used your remaining code since Fourier transform is not right . Because my input signal is a box card function and mathematically i know what its output is going to be that is a sinc function and the way i am using centeredFFT i do get a sync function . > > I need to use FFT and IFFT in my program so that is why i am looking to how to use them correctly hence making a simple test program that i put up where I have a input boxcar and by taking its Fourier Transform I do get a sync function now to see how ifft works i should get back box car function . But i do not . > > the suggestion which you gave does not convince me since Pdft is not a sync function and x axis give me no information about the frequency . > > Could you elaborate your answer ? > > and suggest me possibly How could I recover back Box car after doing ifft with the sync function > > Regards > > > > "Wayne King" <wmkingty(a)gmail.com> wrote in message <hnt0ob$gvm$1(a)fred.mathworks.com>... > > "Saurabh " <saurabhs(a)student.ethz.ch> wrote in message <hnsvof$2f0$1(a)fred.mathworks.com>... > > > Hello > > > > > > I have a rectangular function defined as follows > > > tstop = 10E-3; > > > t=-0.1:dt:0.1; > > > p_i=(1/2)*(sign(t)-sign(t-tstop)); > > > input = figure; > > > plot(t,p_i) > > > rge = 0.05; > > > axis([-rge rge 0 1]); > > > xlabel('time') > > > ylabel('amplitude') > > > title('Input step function') > > > > > > Then I take a Centered FFT > > > > > > [pw,frequencyRange] = centeredFFT(p_i,1/dt); > > > positiveFFT = figure; > > > stem(frequencyRange,abs(pw)); > > > % set(positiveFFT,'Position',[500,500,500,300]) > > > xlabel('Freq (Hz)') > > > ylabel('Amplitude') > > > title('Using the positiveFFT function') > > > grid > > > axis([0 500 0 0.05]); > > > xlabel('frequency') > > > ylabel('amplitude') > > > title('FFT of the input step function') > > > > > > where the function centeredFFT is > > > > > > function [X,freq]=centeredFFT(x,Fs) > > > %this is a custom function that helps in plotting the two-sided spectrum > > > %x is the signal that is to be transformed > > > %Fs is the sampling rate > > > > > > N=length(x); > > > > > > %this part of the code generates that frequency axis > > > if mod(N,2)==0 > > > k=-N/2:N/2-1; % N even > > > else > > > k=-(N-1)/2:(N-1)/2; % N odd > > > end > > > T=N/Fs; > > > freq=k/T; %the frequency axis > > > > > > %takes the fft of the signal, and adjusts the amplitude accordingly > > > X=fft(x)/N; % normalize the data > > > X=fftshift(X); %shifts the fft data so that it is centered > > > > > > Now I want to take IFFT of the function but i do not get back the same input signal > > > > > > orignal_signal = ifft(YfreqDomain); > > > time = ifft(frequencyRange); > > > reterived = figure > > > plot( time,orignal_signal) > > > in > > > > > > Could any one guide in simple terms where I am wrong ? > > > > > > I also tried this > > > > > > orignal_signal = ifft(ifftshift(YfreqDomain)); > > > time = ifft(ifftshift(frequencyRange)); > > > reterived = figure > > > plot( time,orignal_signal) > > > > > > but could not reterive back original signal > > > > Hi, if you want the signal back, why not just > > > > % you never give us dt by the way > > tstop = 10E-3; > > dt=1e-3; > > t=-0.1:dt:0.1; > > p_i=(1/2)*(sign(t)-sign(t-tstop)); > > Pdft = fft(p_i); > > p_m=ifft(Pdft,'symmetric'); > > norm(p_m-p_i) > > > > Wayne Hi, you have to keep in mind that the Fourier transform of a sequence is periodic, so when you say with: Pdft = fft(p_i); plot(abs(Pdft)); You're not seeing sinc-like behavior, that is not correct. You are seeing the sinc-like behavior, you're seeing one period of it. When you use fftshift, you're still seeing one period, you're just looking at the period from a different starting point to ending point. But, it's still the same. Visualize for yourself the periodic extensions--no difference. Wayne Wayne
From: Saurabh on 18 Mar 2010 07:38 Hi Rune Thanks for your reply. But i wonder if i take a some function lets say sine just for case to make myself convinced .. fo = 4; %frequency of the sine wave Fs = 100; %sampling rate Ts = 1/Fs; %sampling time interval t = 0:Ts:1-Ts; %sampling period n = length(t); %number of samples y = 2*sin(2*pi*fo*t); %the sine curve and if you do its Fourier transform using fft YfreqDomain = fft(y); %take the fft of our sin wave, y(t) stem(abs(YfreqDomain)); This doesn’t quite look like what we predicted above. If you notice, there are a couple of things that are missing. The x-axis gives us no information on the frequency. How can we tell that the peaks are in the right place? The amplitude is all the way up to 100 The spectrum is not centered around zero All above goes courtesy to a blog by blinkdagger I am stuck in making some sort of same routine for ifft so that now when I do ifft i could reterive back the starting function ? Regards Rune Allnor <allnor(a)tele.ntnu.no> wrote in message <cc7fb6b1-ff78-460e-bef2-685856aa44b5(a)d27g2000yqf.googlegroups.com>... > On 18 Mar, 11:37, "Saurabh " <saura...(a)student.ethz.ch> wrote: > > > but could not reterive back original signal > > Try this: > > N = 100; > x = randn(N,1); > X = fft(x); > y = real(ifft(X)); > d = max(abs(y-x)) > > d should now contain a very small number (on the order > of 1e-15), which means that y == x, except for numerical > round-off. > > Rune
From: Matt J on 18 Mar 2010 08:12 "Saurabh " <saurabhs(a)student.ethz.ch> wrote in message <hnt3as$qm0$1(a)fred.mathworks.com>... > Hi Rune > > Thanks for your reply. > > But i wonder if i take a some function lets say sine just for case to make myself convinced .. > fo = 4; %frequency of the sine wave > Fs = 100; %sampling rate > Ts = 1/Fs; %sampling time interval > t = 0:Ts:1-Ts; %sampling period > n = length(t); %number of samples > y = 2*sin(2*pi*fo*t); %the sine curve > > > and if you do its Fourier transform using fft > > YfreqDomain = fft(y); %take the fft of our sin wave, y(t) ============== Again, this should be normalized by Fs, YfreqDomain = fft(y)/Fs; %Also, shift it for plotting purposes Also you need to set up your frequency axis, and apply fftshift to center the plot N=n; faxis=( -ceil((N-1)/2):N-1-ceil((N-1)/2) )/N/Ts; stem(faxis,fftshift(abs(YfreqDomain)));
From: Saurabh on 18 Mar 2010 08:47 Hello Matt Thanks for your reply . I run this test code finally to convince myself fo = 4; %frequency of the sine wave Fs = 100; %sampling rate Ts = 1/Fs; %sampling time interval t = 0:Ts:1-Ts; %sampling period n = length(t); %number of samples y = 2*sin(2*pi*fo*t); %the sine curve figure(1) plot(t,y) YfreqDomain = fft(y)./n; %take the fft of our sin wave, y(t) N=n; faxis=( -ceil((N-1)/2):N-1-ceil((N-1)/2) )/N/Ts; figure(2) stem(faxis,fftshift(abs(YfreqDomain))); y_ret = ifft(YfreqDomain).*n; figure(3) plot(t,y_ret) The points I wish to make 1. you suggested to replace N which was length of signal by Fs but i found in the code above no difference. Hence i think normalizing by length of signal makes more meaning 2. Now i can reterive back the signal well It would be instructive if you can point me the difference between DTFT and DFT with some simple example which is implemented in matlab to apreciate the difference Regards "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hnt5ak$sr4$1(a)fred.mathworks.com>... > "Saurabh " <saurabhs(a)student.ethz.ch> wrote in message <hnt3as$qm0$1(a)fred.mathworks.com>... > > Hi Rune > > > > Thanks for your reply. > > > > But i wonder if i take a some function lets say sine just for case to make myself convinced .. > > fo = 4; %frequency of the sine wave > > Fs = 100; %sampling rate > > Ts = 1/Fs; %sampling time interval > > t = 0:Ts:1-Ts; %sampling period > > n = length(t); %number of samples > > y = 2*sin(2*pi*fo*t); %the sine curve > > > > > > and if you do its Fourier transform using fft > > > > YfreqDomain = fft(y); %take the fft of our sin wave, y(t) > ============== > > > Again, this should be normalized by Fs, > > YfreqDomain = fft(y)/Fs; %Also, shift it for plotting purposes > > Also you need to set up your frequency axis, and apply fftshift to center the plot > > N=n; > faxis=( -ceil((N-1)/2):N-1-ceil((N-1)/2) )/N/Ts; > stem(faxis,fftshift(abs(YfreqDomain)));
From: Wayne King on 18 Mar 2010 09:57
"Saurabh " <saurabhs(a)student.ethz.ch> wrote in message <hnt7c9$1od$1(a)fred.mathworks.com>... > Hello Matt > > Thanks for your reply . > > I run this test code finally to convince myself > > fo = 4; %frequency of the sine wave > Fs = 100; %sampling rate > Ts = 1/Fs; %sampling time interval > t = 0:Ts:1-Ts; %sampling period > n = length(t); %number of samples > y = 2*sin(2*pi*fo*t); %the sine curve > figure(1) > plot(t,y) > YfreqDomain = fft(y)./n; %take the fft of our sin wave, y(t) > N=n; > faxis=( -ceil((N-1)/2):N-1-ceil((N-1)/2) )/N/Ts; > figure(2) > stem(faxis,fftshift(abs(YfreqDomain))); > > y_ret = ifft(YfreqDomain).*n; > figure(3) > plot(t,y_ret) > > The points I wish to make > > 1. you suggested to replace N which was length of signal by Fs but i found in the code above no difference. > Hence i think normalizing by length of signal makes more meaning > > 2. Now i can reterive back the signal well It would be instructive if you can point me the difference between DTFT and DFT with some simple example which is implemented in matlab to apreciate the difference > > Regards > > "Matt J " <mattjacREMOVE(a)THISieee.spam> wrote in message <hnt5ak$sr4$1(a)fred.mathworks.com>... > > "Saurabh " <saurabhs(a)student.ethz.ch> wrote in message <hnt3as$qm0$1(a)fred.mathworks.com>... > > > Hi Rune > > > > > > Thanks for your reply. > > > > > > But i wonder if i take a some function lets say sine just for case to make myself convinced .. > > > fo = 4; %frequency of the sine wave > > > Fs = 100; %sampling rate > > > Ts = 1/Fs; %sampling time interval > > > t = 0:Ts:1-Ts; %sampling period > > > n = length(t); %number of samples > > > y = 2*sin(2*pi*fo*t); %the sine curve > > > > > > > > > and if you do its Fourier transform using fft > > > > > > YfreqDomain = fft(y); %take the fft of our sin wave, y(t) > > ============== > > > > > > Again, this should be normalized by Fs, > > > > YfreqDomain = fft(y)/Fs; %Also, shift it for plotting purposes > > > > Also you need to set up your frequency axis, and apply fftshift to center the plot > > > > N=n; > > faxis=( -ceil((N-1)/2):N-1-ceil((N-1)/2) )/N/Ts; > > stem(faxis,fftshift(abs(YfreqDomain))); Hi Saurabh, You don't need to use fftshift() if you don't want to. fo = 4; %frequency of the sine wave Fs = 100; %sampling rate Ts = 1/Fs; %sampling time interval t = 0:Ts:1-Ts; %sampling period n = length(t); %number of samples y = 2*sin(2*pi*fo*t); %the sine curve freq = 0:(Fs/length(y)):Fs/2; YfreqDomain = fft(y); stem(freq,2/length(y)*abs(YfreqDomain(1:length(y)/2+1))) figure; y_ret=ifft(YfreqDomain,'symmetric'); plot(t,y_ret); In MATLAB you are always computing the DFT, the DTFT is a continuous function of frequency. You are always evaluating the Fourier transform of a sequence on a discrete grid. Wayne |