Monday, October 19, 2009
Thursday, October 15, 2009
FFT in Matlab, from WIKI
fft - Discrete Fourier transform
Syntax
Y = fft(X)
Y = fft(X,n)
Y = fft(X,[],dim)
Y = fft(X,n,dim)
Definition
The functions Y=fft(x) and y=ifft(X) implement the transform and inverse transform pair given for vectors of length 
where
is an 
Description
Y = fft(X) returns the discrete Fourier transform (DFT) of vector X, computed with a fast Fourier transform (FFT) algorithm.
If X is a matrix, fft returns the Fourier transform of each column of the matrix.
If X is a multidimensional array, fft operates on the first nonsingleton dimension.
Y = fft(X,n) returns the n-point DFT. If the length of X is less than n, X is padded with trailing zeros to length n. If the length of X is greater than n, the sequence X is truncated. When X is a matrix, the length of the columns are adjusted in the same manner.
Y = fft(X,[],dim) and Y = fft(X,n,dim) applies the FFT operation across the dimension dim.
Examples
A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal. Consider data sampled at 1000 Hz. Form a signal containing a 50 Hz sinusoid of amplitude 0.7 and 120 Hz sinusoid of amplitude 1 and corrupt it with some zero-mean random noise:
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)') 
It is difficult to identify the frequency components by looking at the original signal. Converting to the frequency domain, the discrete Fourier transform of the noisy signal y is found by taking the fast Fourier transform (FFT):
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)|') 
The main reason the amplitudes are not exactly at 0.7 and 1 is because of the noise. Several executions of this code (including recomputation of y) will produce different approximations to 0.7 and 1. The other reason is that you have a finite length signal. Increasing L from 1000 to 10000 in the example above will produce much better approximations on average.
Algorithm
The FFT functions (fft, fft2, fftn, ifft, ifft2, ifftn) are based on a library called FFTW [3],[4]. To compute an 
When 
For most 
The execution time for fft depends on the length of the transform. It is fastest for powers of two. It is almost as fast for lengths that have only small prime factors. It is typically several times slower for lengths that are prime or which have large prime factors.
Friday, September 25, 2009
g++ extra qualification
很多比较旧的代码会有如写法
class Foo
{
int Foo::Foo(void);
}
在g++ 4.1以后会报错 extra qualification
直接改为
class Foo
{
int Foo(void);
}
即可
The author`s blog: http://hi.baidu.com/zjugator/blog/item/77bb6cec02bee22163d09f7d.html
Sunday, September 13, 2009
GCC configureation
gcc App.o -o App -L/usr/local/lib -L/usr/lib -L/usr/lib64 -DNO_SDL_MAIN -lz -lzip -lzipd -ldl -lX11 -lXext -lpthread -lGLU -lGL -ljpeg -lpng -lavformat -lavcodec -lavutil -lavdevice -lSDL -lG3D -lGLG3D -lG3Dd -lGLG3Dd -lm
这是GCC编译G3D例子starter的输入！看上去确实吓人，在刚接触GCC的时候，这个东西几乎让我疯狂了！在抓狂了一周后，今天凌晨终于解决了编译和链接的问题！现在谈谈这两天的感觉。希望抛砖引玉！
1. Linux下的library 和header都是搜索式的，比如有两个路径，/usr/local/lib和/lib，里面分别有一个叫G3D的header，当你使用linux下的编译器时，比如gcc，他会先搜索一个路径，在搜索另一个路径，先到先的，这个顺序是取决于你的输入顺序或设置的，一些新安装的库可能需要你对输入做一定的排序才能正确使用。G3D的安装过程中，我就犯了一个这样的错误，/usr/local/lib首先搜索，其下的libzipd.a没有更新，所以不管我怎么修改配置，重新编译，更换版本都无济于事，最后的解决办法很简单，要么删除/usr/local/lib下的库或者将新编译的库文件设置到较高的搜索优先级。Windows下就相对简单了，他们是登记式的，或者是更细化的搜索，你只要在注册表，环境变量或编译器中登记要使用的类库，虽然这个过程linux也有，但windows是没有固定的/lib和/include目录的，就是没有默认搜索，所以windows这个过程更直观一些。
2.“Undefined reference to .......” in link step. 这个问题通常因为library没有找到，比如你安装library的路径没有被编译器搜索，而include却被搜索到了，所以就出现了编译成功而链接出错！这里更需要注意的是gcc与window的错误信息是有区别的，在windows上，undefind refernce 肯定出现在编译阶段，链接阶段会说undefind symbol...如果你习惯了vc的错误信息，很容易把gcc link的错误当作编译的错误，所以不管你怎么改，都是不会有效的......
3.gcc的库是搜索式的，所以你不需要定义库的位置，只需要定义库的名字，如果你的安装路径不是默认的，你还需要设定一个搜索路径，gcc会一个接一个的搜索，当然，这里也会有谁先谁后的问题！library的命名方式是lib开头，然后是库的名字，最后的后缀是lib。gcc输入这是-l开头，紧跟着是lib的名字。比如，我向gcc传入 -lG3D， 在../lib中，他应该叫做libG3D.lib。最后用-lm表示library输入结束，而且要在最后！-L说明搜索的路径，FIFS！
这是GCC编译G3D例子starter的输入！看上去确实吓人，在刚接触GCC的时候，这个东西几乎让我疯狂了！在抓狂了一周后，今天凌晨终于解决了编译和链接的问题！现在谈谈这两天的感觉。希望抛砖引玉！
1. Linux下的library 和header都是搜索式的，比如有两个路径，/usr/local/lib和/lib，里面分别有一个叫G3D的header，当你使用linux下的编译器时，比如gcc，他会先搜索一个路径，在搜索另一个路径，先到先的，这个顺序是取决于你的输入顺序或设置的，一些新安装的库可能需要你对输入做一定的排序才能正确使用。G3D的安装过程中，我就犯了一个这样的错误，/usr/local/lib首先搜索，其下的libzipd.a没有更新，所以不管我怎么修改配置，重新编译，更换版本都无济于事，最后的解决办法很简单，要么删除/usr/local/lib下的库或者将新编译的库文件设置到较高的搜索优先级。Windows下就相对简单了，他们是登记式的，或者是更细化的搜索，你只要在注册表，环境变量或编译器中登记要使用的类库，虽然这个过程linux也有，但windows是没有固定的/lib和/include目录的，就是没有默认搜索，所以windows这个过程更直观一些。
2.“Undefined reference to .......” in link step. 这个问题通常因为library没有找到，比如你安装library的路径没有被编译器搜索，而include却被搜索到了，所以就出现了编译成功而链接出错！这里更需要注意的是gcc与window的错误信息是有区别的，在windows上，undefind refernce 肯定出现在编译阶段，链接阶段会说undefind symbol...如果你习惯了vc的错误信息，很容易把gcc link的错误当作编译的错误，所以不管你怎么改，都是不会有效的......
3.gcc的库是搜索式的，所以你不需要定义库的位置，只需要定义库的名字，如果你的安装路径不是默认的，你还需要设定一个搜索路径，gcc会一个接一个的搜索，当然，这里也会有谁先谁后的问题！library的命名方式是lib开头，然后是库的名字，最后的后缀是lib。gcc输入这是-l开头，紧跟着是lib的名字。比如，我向gcc传入 -lG3D， 在../lib中，他应该叫做libG3D.lib。最后用-lm表示library输入结束，而且要在最后！-L说明搜索的路径，FIFS！
Tuesday, September 1, 2009
Least Squares Approximation
| Suppose you are in a science class and you receive these instructions:
Find the temperature of the water (in degrees Celsius) at the times 1, 2, 3, 4, and 5 seconds after you have applied heat to the container. Conduct your experiment carefully. Graph each data point with time on the x-axis and temperature on the y-axis. Your data should follow a straight line. Find the equation of this line. The data from the experiment looks like this when charted and graphed:
This is a common problem with experimental sciences because the data points that we measure seldom fall on a straight line. Therefore, scientists try to find an approximation. In this case, they would try to find the line that best fits the data in some sense. The first problem is to define "best fit." It is convenient to define an error as a distance from the actual value of y for x (the value that was measured in the experiment) to the predicted value of y for x. Therefore, it seems reasonable that the "best fit" line would somehow minimize the errors, but how? You could minimize the sum of the absolute values of the errors; this is called the You have just read a lot of new information, so let's illustrate the concepts with our example. We have the graph of the data above. Now we need to guess which line best fits our data. If we assume that the first two points are correct and choose the line that goes through them, we get the line y = 1 + x. If we substitute our points into this equation, we get the following chart. The points and line are graphed below.
If we choose the line that goes through the points when x = 3 and 4, we get the line y = 4 + x. Will we get a better fit? Let's look at it.  The sum of the squares of the error is 18. That is a better fit, but can we do even better? Let's try the line that is half way between these two lines. The equation would be y = 2.5 + x. It looks like this:
A line in slope-intercept form looks like c0 + c1x = y where c0 is the y -intercept and c1 is the slope. We want to find c0 and c1 such that c0 + c1xi = yi is true for all our data points:
 In general, we cannot solve this system because the system is usually inconsistent because it is overdetermined. In other words, we have more equations than unknowns (the unknowns are the two variables, c0 and c1, for which we are trying to solve). There is a system of equations called the normal equations that can be used to find least squares solution to systems with more equations than unknowns.
Theorem 8.1 Let X be an m by n matrix such that XTX is invertible, then the soltution to the normal equations, XTXc = XTy, is the least squares approximation to c in Xc = y.
Remark 25 It is important to remember that the solution to the normal equations is only an approximation to c for Xc = y. It is not equal to c because Xc = y is inconsistent, so it has no solution. In other words, there does not exist a vector, c, that makes Xc = y a true statement. Therefore, we use the normal equations to approximate c.
Remark 26 For now, you don't need to check to see if XTX is invertible because most of the systems that we encounter will meet this requirement. However, if you cannot find a solution to the normal equations, you should check to see if XTX is invertible. The normal equations will give us the "best fit" line (or curve) every time according to the way we defined "best fit." The proof of this is at the end of this chapter. Let's try applying the normal equations to our system. First, we multiply so that we have a system that we can solve.
 The sum of the squares of the error is 2.7. This is a great improvement over our guesses and we know that we cannot do any better. In general, if we have n data points, we solve XTXc = XTy with
 ,  and   or  tell you to continue in the same pattern. What if we are told that our data is not supposed to fit a straight line, but instead falls in the shape of a parabola? Consider the following data from another experiment:
 We can find the curve that best fits our data in a similar manner. The general equation for a parabola is c0 + c1x + c2x2 = y. Therefore, we want to find the values of the coefficients, c1 c2, and c3, so that the curve we find best fits these equations:
 These coefficients indicate that the curve we want is  Let's graph this curve and fill in our chart:
 We find that the sum of the squared errors is
 Notice that the normal equations used to find the best fit line and the best fit parabola have the same form. Do you think that we could expand this to higher degree polynomials? Yes, we can. In general, we use the normal equations XTXc = XTy with
 where m represents the degree of the polynomial curve that you wish to fit and n represents the number of data points. The least squares "best fit" curve for these equations is c0 + c1x +c2x2 + … + cm - 1xm - 1 + cmxm. Remember that the degree is the highest power of the variable in your equation. A line is a first degree polynomial and a parabola is a second degree polynomial. If we can find the best fit curve for any degree polynomial, why don't we always use a higher degree polynomial and fit the data better? After all, if we have n data points and fit them to a polynomial of degree n - 1, we will have a perfect fit every time because our systems would not be inconsistent. However, our goal is not just to find a curve that fits the data closely. Usually, we want the curve to predict what would happen between our data points. If we choose a curve that exactly fits all our data points, we are incorporating the error in our measurements into our model unless the model fits the data exactly (which occurs only rarely). Unfortunately, there is no set rule for deciding what degree polynomial should be used to fit the data. However, first and second degree polynomials provide the simplest models and should fit most of your data until you start modeling more complicated systems. If you notice, we said that we usually fit a curve so that we can predict what would happen between our data points. Predicting an outcome between data points is called interpolation. Why didn't we say anything about predicting the behavior beyond our data points? Predicting an outcome beyond the data is called extrapolation. It is usually dangerous to extrapolate much beyond the data because we have no indication that the data will continue to follow the same curve since our curve was only fit to the data. For example, we measured the height of a teenage boy every year for a few years and charted his growth. The growth appeared linear, so we fit a line to the data and got y = 32 + 2.25x. We have graphed the data with age on the x-axis and height on the y-axis.
 If we extrapolate back several years, this young man was over two and a half feet tall when he was born. According to this model, he will never stop growing, so he will be 8 feet 4 inches tall by the time he is 30 and almost 14 feet tall by the time he is 60. Do you think that this is an accurate prediction? If the temperature at the airport on the 4th of July was in the 90's for two years in a row, would it be reasonable to predict that the temperature in January between those years was also in the 90's? No, it would not. We have two problems with this model. One problem is that we only have 2 data points. You can always find a line that fits the two points, but there is no reason to believe that the relationship between the day of the year and the temperature is a linear relationship. Also, we didn't take into account other factors that could affect our model such as the pattern of the seasons. These are problems that can arise when you model a situation. When we start modeling situations and using least squares to make predictions, we are entering the world of statistics. That means that we must think about what the data represents rather than just apply the normal equations. There are many interesting applications of statistics that you can explore in another course. However, using matrices, you already know one way to find a "best fit" curve for your data. |
fit. It would also be reasonable to find the biggest error for each line and choose the line that minimizes this quantity; this is called the 
fit. However, the fit that is used most often is the 
fit which is called the 
and 
Using our definition of least squares "best fit," you will not be able to find a parabola that fits the data better than this one. In general, to find the parabola that best fits the data, you use the normal equations XTXc = XTy with
Subscribe to:
Posts (Atom)