Author |
Topic |
|
ubitelli@n1
Brazil
8 Posts |
Posted - 05/19/2003 : 4:03:11 PM
|
Hi people.
Regard the last post, I did not find the book Hans suggested but I have foud others books. However I still have questions about the FFT smooth procedure. If I am not wrong, the procedure is something like this: 1)Take the FFT of the data. 2)Multiply this FFT by a function (a filter?) that has a value of 1 at zero Hz and 0 at the frequency of 1/(N*dx) where N is the number of points used to smoothing and dx is the interval in the x axis. In this way, we remove the components above the frequency cutoff of 1/(N*dx). 3)Take the inverse FFT.
Well, in the item 2, are the data above cutoff removed or replaced by zero? I think they must be replaced because the number of data points must not change.
Anyway, I used this procedure with simulated data and my results did not agree with Origin (I did´t use any function but replaced the data by zero using brute force - with few points it is easy).
What is wrong? Thanks Ricardo.
|
|
easwar
USA
1964 Posts |
|
ubitelli@n1
Brazil
8 Posts |
Posted - 05/20/2003 : 1:56:06 PM
|
Hello Easwar
You are right. The cutoff can´t be so sharp. Yes Easwar, I would like to know how this function is constructed. Could you send me some information on this? Thanks Ricardo. |
|
|
easwar
USA
1964 Posts |
Posted - 05/20/2003 : 4:56:52 PM
|
Hi Ricardo,
The cutoff frequency is determined using the formula: FCutoff = 1 / (msmooth * xinc) where msmooth is the number of smooth points you specify, and xinc is the increment/step value in time. xinc is computed as the difference between the first two x values.
A parabola is then generated that is symmetric about 0, and has a value of 1 at freq=0 and has a value of 0 at freq=FCutoff.
This parabola is used to filter/multiply the FFT spectrum, and the inverse FFT computation is performed to obtain the smoothed data.
If you are trying to code this yourself, and you want to play with the parabola/filter shape etc, I would recommend using the NAG fft functions - you can use the FFT() function under matrix.h - I am assuming here that you have version 7.
We will add this information to our documentation.
Easwar OriginLab.
|
|
|
ubitelli@n1
Brazil
8 Posts |
Posted - 05/22/2003 : 12:33:49 PM
|
Hi Easwar
I have made some calculations using your hints but the things aren´t still good. Is the parabolic function just a function of this type: 1-(1/Fc^2)*x^2 where Fc is the cutoff frequency and x are the x axis values? If this function is correct something is wrong with my procedure. What I am doing is taking the FFT of the y data, multiplying this FFT by the function defined above and taking the inverse FFT. The results don´t match with the Origin´s results. If you want I can send you the file with the data points. I am using only 128 points.
Thanks Ricardo.
PS.: I have Origin version 6.1
|
|
|
ubitelli@n1
Brazil
8 Posts |
Posted - 05/30/2003 : 2:06:24 PM
|
Any body here? Nothing new on the FFT smoothing procedure? What a hell of parabolic function!!!
I will try suicide.....
Ricardo. |
|
|
easwar
USA
1964 Posts |
Posted - 05/30/2003 : 4:54:12 PM
|
Hi Ricardo,
Sorry for the delay in replying. I presume you are trying to do this using the GUI by doing FFT, then multiplying with filter etc. You need to be careful in doing this because you need to multiply both positive and negative frequencies by the parabolic shape.
Try the LabTalk code posted below. This code assumes that your time data is in data1_a, and the signal data is in data1_b. Running the script will then create FFT1 and IFFT1 worksheets. It will also bring up a dialog asking you for number of smoothing points. The smoothed curve will be in IFFT1_Real.
Note that the result from this script is slightly different from what the built-in tool does. This is because the tool internally pads the data to 2^n number of points and then does FFT etc. But this script could be a good starting point if you want to change/improve the filter.
One way to run this script is to edit custom.ogs file on your root Origin directory, then comment out lines below [Main] section, and paste this code under there. Then you can run code by clicking on the Custom Routine button on the Standard Toolbar.
In version 7, Origin is shipped with a selection of NAG numerical libraries, and also supports C programming - so filter code can be written using C and one can call NAG library for performing FFT etc.
Hope this helps you.
Easwar OriginLab.
// First check if worksheets FFT1 and IFFT1 exist - if yes, clear them if(exist(FFT1) == 2) clearworksheet FFT1; if(exist(IFFT1) == 2) clearworksheet IFFT1;
// Perform forward FFT fft.reset(); fft.forward = 1; fft.forward.timeData$ = Data1_A; fft.forward.tdelta = Data1_A[2] - Data1_A[1]; fft.forward.realData$ = Data1_B; window -t W fft FFT1; fft.output.samplingdata$ = FFT1_Freq; fft.output.realdata$ = FFT1_Real; fft.output.imagdata$ = FFT1_Imag; fft.output.ampdata$ = FFT1_r; fft.output.phasedata$ = FFT1_Phi; fft.output.powerdata$ = FFT1_Power; fft.real = 1; fft.normalize = 0; fft.shifted = 0; fft.windowing = 1; fft.spectrum = 1; fft.unwrap = 0; fft.forward();
// Bring up dialog and ask user for number of points for smoothing getn "Enter number of smooth points;" nPts; get Data1_A -e ii; // Quit if number of points don't make sense if( (nPts < 1) || (nPts > (ii/2)) ) { type Num. of points incorrect!; break; }
// Compute cutoff frequency FCutoff = 1 / ( nPts * (Data1_A[2] - Data1_A[1]) ); // Find index in FFT dataset where cutoff frequency occurs ii1 = xindex(FCutoff, FFT1_Real); // Find length of FFT dataset and find mirror point get FFT1_Freq -e ii; ii2 = ii - ii1;
// Create a Filter dataset create Filter ii; set Filter -e ii; // Fill it first with all 1s Filter = 1;
// Loop thru filter and create parabolic shape loop(jj, 2, ii1) { Filter[jj] = sqrt(1 - jj^2/ii1^2); Filter[ii-jj+2] = sqrt(1 - jj^2/ii1^2); }
// Set middle values in filter to zero loop (jj, ii1, ii2) { Filter[jj] = 0; }
// Multiply FFT real and imaginary with Filter shape FFT1_Real *= Filter; FFT1_Imag *= Filter;
// Delete Filter dataset delete Filter;
// Perform backward FFT fft.reset(); fft.forward = 0; fft.backward.freqData$ = FFT1_Freq; fft.backward.fDelta = FFT1_Freq[2] - FFT1_Freq[1]; fft.backward.realData$ = FFT1_Real; fft.backward.imagData$ = FFT1_Imag; window -t W fft IFFT1; fft.output.samplingdata$ = IFFT1_Freq; fft.output.realdata$ = IFFT1_Real; fft.output.imagdata$ = IFFT1_Imag; fft.output.ampdata$ = IFFT1_r; fft.output.phasedata$ = IFFT1_Phi; fft.output.powerdata$ = IFFT1_Power; fft.real = 1; fft.normalize = 0; fft.shifted = 0; fft.windowing = 1; fft.spectrum = 1; fft.unwrap = 0; fft.backward(); // End of script file
|
|
|
ubitelli@n1
Brazil
8 Posts |
Posted - 06/03/2003 : 11:41:32 AM
|
Hi Easwar
Thanks for your code. I have runned in Origin 6.1 and the results are slightly different as you said. But now I know how the calculations are done. The "trick" was the negative frequency part of the spectrum isn´t it?
Thank you. Ricardo. |
|
|
rtoomey
USA
184 Posts |
|
|
Topic |
|