Guy and Seth on Simulink
January 30th, 2009
MathWorks Conversations and the FFT
If you are like me, you read the doc, a lot. I am often
clicking on the help just to verify my understanding of a function’s syntax, or
the behavior of a block. That is why I wasn’t surprised to find myself
involved in a detailed discussion about the documentation for FFT.
You go it… just another day at the MathWorks. For some
context, the doc example generates a signal corrupted with noise, and then uses
the FFT to extract the frequency components.
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t)); % Sinusoids plus noise
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

There were two questions that started the discussion. The
first one,
Why does the FFT example in the doc have so much code?
We can ignore that for now and focus on the second question,
and the real meat of the discussion:
Why does the FFT example result in an amplitude of 1?
My colleague Jeoffrey surprised me with his intimate
familiarity with the discrete time Fourier transform and his ability to produce
the complicated equations out of thin air. I told him, “I want that in a blog
post!” So here it is. Today’s featured guest blogger is Jeoffrey Young:

Why does the FFT example result in an amplitude of 1? By Jeoffrey Young
In this post I will talk about one way of looking at where
the scaling comes from in
the following example from Matlab 7.7 (R2008b)’s help documentation:

To start with, let’s remember that FFT is simply the sampled
Discrete Time Fourier Transform of a signal. We know from DTFT that (I’ll use
the cosine function here for simplicity):

where and (Ever wonder
why this looks exactly the same as the Fourier Transform of the continuous time
cosine signal? More on this later.) However, it is not enough to sample a
continuous time signal that is not time-limited, such as a cosine function.
Truncating so that we
can represent the signal on a computer is equivalent to multiplying the time-domain
signal by the following rectangular window:

The DTFT of this signal is the Digital Sinc Function:

Note that by
L’Hopital’s Rule.
Thus, from the time-domain multiplication property of
Fourier Transform:

Here we see that the DTFT of an L-point truncated cosine function
has a maximum amplitude of (if we ignore
the effects of aliasing from the overlap and the fact that DTFT is periodic).
Question about the scale factor used with FFT are also sent
to technical support. This happens often enough that tech support posted a
solution titled “Why
is the example of using an FFT to get a frequency spectrum scaled the way it
is?” We even have a technical note on the general topic of spectral
analysis with FFT: Tech Note:
1702 - Using FFT to Obtain Simple Spectral Analysis Plots. As we see it used in this post, the section below describes how the scale factor may also depend on the signal of interest.
Why does looks
exactly like ?
Remember the following relationship between the DTFT of the
sampled analog signal and the FT of the original signal:

In other words, the Discrete Time Fourier transform is
scaled by . However, the
DTFT of a sinusoid is the Delta Function, which is defined more by its
properties than its value. Thus, its amplitude has to be scaled when the axis
is modified :

So that:

Now it's your turn
If you have any
thoughts on how scaling the FFT results may be useful, please post your comments here.
12:27 UTC |
Posted in Fundamentals, Guest Blogger, Signal Processing |
Permalink |
10 Comments »
You can follow any responses to this entry through the RSS 2.0 feed.
You can skip to the end and leave a response. Pinging is currently not allowed.
Leave a Reply
|
Sir,
I have tried a lot but i am unable to find fft in simulunk. Can you tell me which block should i use to find fft and which scope block to see it?
Thankyou.
@Saurabh – The blocks I think you are looking for are part of the Signal Processing Blockset. You can find an FFT block in the Power Spectrum Estimation library. There is also a Spectrum Scope in the Signal Processing sinks. If you don’t have the Signal Processing Blockset, output your signal data to the MATLAB work space and compute and FFT from there. I hope that helps.
How can I have same example but instead AC(1 to 10V 50 or 60Hz) and DC(0.5 to 10 V) then adding AC+DC but out put should be only DC signal? I will be very glad to know if this kind of*. m code can be transferred to simulink block where the block accepts AC+DC as input and the out put as DC signal only.
This example is quite good to know.
Thanks for this. It answers most of my questions. Can you recommend a good, IN PRINT reference text that covers all of these issues? Another one that isn’t well covered is how to convert DFT’s between “normalized” and un-normalized frequencies. The reference texts you give are all dated and mostly out of print; or they cover one issue but not the others.
Thanks again.
@Jay – I’m glad this helped! Any Digital Signal Processing text should cover this material. For example, the text “Discrete-Time Signal Processing by Oppenheim, A. V. and R. W. Schafer” should still be in print. As for the different time a frequency terminologies used in MATLAB, the following help documentation provides the complete list:
http://www.mathworks.com/access/helpdesk/help/toolbox/dspblks/ug/f13-54010.html#f13-23126
Dear Seth,
I’m starting with just a frequency response and I’m having trouble re-adjusting the response( ie scaling) such that ifft could be used to generate a time domaine response. Here is a sample code that I’ve generated for a one simple transfer function example(ie g(jw)=1/(jw+1)).
k=1; p=-1; dt=1e-2; tstop=30.01; %Time period t=0:dt:tstop; %Apply Freq Response freq=0:1/tstop:1/dt; jw=j*2*pi*freq; w=2*pi*freq; for d=1:length(freq) g(d)=k/(jw(d)-p); end % %double sided response n=length(g) m=fliplr(conj(g(2:((n)/2+1)))); mm=[g(1:(n)/2) m]*n; x=real(ifft(mm))/tstop; plot(t,x) hold on %Known Time domain result from Laplace for d=1:length(t) h(d)=exp(p*t(d)); end plot(t,h,'r')How would I handle the maximum frequency point with an even number of points? Is my scaling of my result correct? I hope that you can help me clear the issues behind ifft/fft. Thanks once again,
Dragan
Hi Seth,
I am having a scaling problem with performing FFT in Simulink as well. When the output is displayed on a vector scope the resulting magnitude is too high compared to the theoretical value. While this alone is not the problem. The amplitude varies with FFT length and with the stop time too.
Any suggestions on scaling it appropriately?
Regards,
Chai
Hi Seth,
Thanks for the naswer regarding Saurabh’s question. I’am having the same problem. I don’t have the Signal Processing Blockset install in my Matlab. I wanted to write a script to generate it form my output signal. As a newbie with Matlab I don’t really know how to introce the time and the different part of the script. This is what I have writen but it doesn’t work.
dt = 1/100;
T = 100;
t = (0:dt:T-dt);
x1 = voltage2.signals.values(:,1);
y = fft(x1)/length(x1)*2;
m = abs(y);
m(1)= m(1)/2;
f = (0:length(y)-1)*100/length(y);
plot(f,20*log10(m)),
hold on
ylabel(‘Abs. Magnitude dB’),
x2 = voltage3.signals.values(:,1);
y = fft(x2)/length(x2)*2;
m = abs(y);
m(1)= m(1)/2;
f = (0:length(y)-1)*100/length(y);
plot(f,20*log10(m),’-.’),
axis ([0 max(f)/2 -60 10])
hold off
h = legend(‘Uu’,'Uu0′);
It’s the ouput voltage of a 2 level 4 quadrant IGBT circuit.
X1 is the line to virtual neutral of the “Y” connected out put load voltage and X2 is the line to ground load voltage.
Please could u help me.
Many thanks.
Hello,
I think it would be useful to have a matlab function that directly computes time domain response from s-domain system transfer function available only as numerical data:
G(s=jw). While I realize there are subtleties in the FFT and iFFT scaling I’d rather not be bothered by these issues, since my focus is on dyanmics rather than signal procressing techniques and algorithms.
Amit.
Dear Sir,
I have a problem at hand that requires urgent attention. Thing is I’ve never used Simulink before and right now its required for me to FFT a continuous signal output from my vehicle model.
Thing is it is a continuous signal and I am facing problems trying to even implement the FFT. Even when there is no errors there would not be any output. So yes i’m desperate for help.
thanks.
XY