2.13.4 Signal ProcessingSignal-Processing
Origin provides a collection of X-functions and LabTalk functions for signal processing, ranging from smoothing noisy data to Fourier Transform (FFT), Inverse Fourier Transform (IFFT), Short-time FFT, Convolution and Correlation, FFT Filtering, and Wavelet analysis.
The X-Functions are available under the Signal Processing category and can be listed by typing the following command:
lx cat:="signal processing*";
Some functionality such as Short-time FFT and Wavelets are only available in OriginPro.
The LabTalk functions for signal processing are available under Signal Processing Functions, which are used to separately compute FFT results such as amplitude, phase, etc., for multiple datasets and arrange their results in desired columns. Meanwhile, one can utilize these functions to compare DC offset removed magnitudes among different datasets.
The following sections provide some short examples of calling the signal processing X-Functions and LabTalk functions from script.
Smoothing
Smoothing noisy data can be performed by using the smooth X-Function.
// Smooth the XY data in columns 1,2 of the worksheet
// using SavitzkyGolay method with a third order polynomial
range r=(1,2); // assume worksheet active with XY data
smooth iy:=r meth:=sg poly:=3;
To smooth all plots in a layer, you can loop over the plots as below:
// Count the number of data plots in the layer and save result in
//variable "count"
layer -c;
// Get the name of this Graph page
string gname$ = %H;
// Create a new book named smooth - actual name is stored in bkname$
newbook na:=Smoothed;
// Start with no columns
wks.ncols=0;
loop(ii,1,count) {
// Input Range refers to 'ii'th plot
range riy = [gname$]!$(ii);
// Output Range refers to two, new columns
range roy = [bkname$]!($(ii*2-1),$(ii*2));
// Savitsky-Golay smoothing using third order polynomial
smooth iy:=riy meth:=sg poly:=3 oy:=roy;
}
FFT and IFFT
The following example shows how to individually obtain the FFT results by LabTalk functions listed under Signal Processing Functions.
newbook;
// Import the signal data
string fname$ = system.path.program$ + "Samples\Signal Processing\fftfilter1.DAT";
impASC fname:=fname$;
// Add 5 columns to store different quantities
worksheet -a 5;
// Set column label to distinguish
col(B)[L]$ = "Raw Signal";
// Set data type to complex prior to store complex results from FFT analysis
wks.col3.numerictype = 11;
col(C)[L]$ = "FFT Complex";
// Calculate FFT complex results
col(C) = fftc(col(B)); //with no shift
col(D)[L]$ = "Frequency";
wks.col4.type = 4; //Set X column type
// Compute the frequencies based on time data in column A
col(D) = fftfreq(col(A)[2]-col(A)[1], wks.col1.nrows);
col(E)[L]$ = "Magnitude";
// Get FFT magnitude results
col(E) = fftmag(col(C)); // Use fftmag(col(C), 2) to obtain Two-Sided and Shifted magnitude
col(F)[L]$ = "Phase";
col(F) = fftphase(col(C)); // Use fftphase(col(C), 2, 1, 0) to obtain Two-sided, unwrapped phase with radian unit
wks.col7.numerictype = 11;
col(G)[L]$ = "Shifted FFT Complex";
col(G) = fftshift(col(C));
// Update sparklines for calculated results
sparklines sel:=0 c1:=4 c2:=6;
// Auto adjust column width to fit content
wautosize;
The following example shows how to obtain the IFFT result from shifted FFT complex result.
//Continue from the example above, suppose we have shifted FFT results in column G
// Add 2 more columns to store different quantities
worksheet -a 2;
// Label columns to distinguish
col(H)[L]$ = "Unshifted FFT Complex";
// Set the data type to be complex
wks.col8.numerictype = 11;
// Unshift the shifted FFT results before doing IFFT
col(H) = ifftshift(col(G));
col(I)[L]$ = "IFFT result";
wks.col9.numerictype = 11;
// Compute inverse FFT from unshifted FFT results
col(I) = invfft( col(H) );
The following example shows how to specify a window in FFT analysis.
newbook;
// Import the signal data
string fname$ = system.path.program$ + "Samples\Signal Processing\Chirp Signal.dat.";
impASC fname:=fname$;
// Define a range variable for input signal
range rSignal = col(B);
int wSize = rSignal.nrows;
col(C)[L]$ = "Apply Blackman window";
col(C) = col(B) * windata(6, wSize); // Apply Blackman window
col(D)[L]$ = "Amplitude";
col(D) = fftamp( fftc(col(C)) ); //Without Origin window correction
The following example shows how to perform fft analysis on multiple columns using loop and arrange obtained amplitude and phase columns side by side. You can also calculate FFT results for multiple columns directly on worksheet by selecting multiple columns on the worksheet, and right click to select Set Multiple Columns Values....
//Prepare multiple column data for FFT analysis
newbook;
int nc = 5, ii;
wks.ncols = 4*nc+2;
//Fill x column with time data
col(A) = data(0, 1-1e-3, 1e-3);
//Fill five columns with sum of an f1 Hz sinusoid and an f2 Hz sinusoid
for( ii = 1; ii<=nc; ii++ )
{
double f1, f2;
f1 = 50 + 20*ii;
f2 = 100 + 20*ii;
wcol(ii+1)[L]$ = "Data$(ii)";
wcol(ii+1) = (2+ii)*sin(2*pi*f1*col(A))+(8-ii)*sin(2*pi*f2*col(A))+rnd();
}
//Subtract mean before FFT to remove DC offset
for( ii = 1; ii<=nc; ii++ )
{
wcol(nc+ii+1)[C]$ = "Subtract mean of Data$(ii)";
wcol(nc+ii+1) = wcol(ii+1)-mean( wcol(ii+1) );
}
//Calculate frequency for FFT results
wks.col$(2*nc+2).type = 4; //Set X column type
wcol(2*nc+2)[L]$ = "Frequency";
wcol(2*nc+2) = fftfreq( col(A)[2]-col(A)[1], wks.col1.nrows );
//Calculate FFT amplitude and phase for five columns
for( ii = 1; ii<=nc; ii++ )
{
//FFT amplitude result
wcol(2*nc+ii+2)[L]$ = "Amplitude";
wcol(2*nc+ii+2)[C]$ = "Data$(ii)";
wcol(2*nc+ii+2) = fftamp( fftc( wcol(nc+ii+1) ) );
//FFT phase result
wcol(3*nc+ii+2)[L]$ = "Phase";
wcol(3*nc+ii+2)[C]$ = "Data$(ii)";
wcol(3*nc+ii+2) = fftphase( fftc( wcol(nc+ii+1) ) );
}
//Plot FFT amplitude and phase results for five columns in two layers
plotstack iy:=($(2*nc+2),$(2*nc+3):$(4*nc+2)) portrait:=0 order:=0 layer:=2 number:="5 5";
FFT and Filtering
The following example shows how to perform 1D FFTFast Fourier Transform of data using the fft1 X-Function.
// Import a sample file
newbook;
fname$ = system.path.program$ + "Samples\Signal Processing\fftfilter1.dat";
impasc;
// Perform FFT and get output into a named tree
Tree myfft;
fft1 ix:=2 rd:=myfft rt:=<none>;
// You can list all trees using the command: list vt
Once you have results in a tree, you can do further analysis on the output such as:
// Copy desired tree vector nodes to datasets
// Locate the peak and mean frequency components
dataset tmp_x=myfft.fft.freq;
dataset tmp_y=myfft.fft.amp;
// Perform stats and output results
percentile = {0:10:100};
diststats iy:=(tmp_x, tmp_y) percent:=percentile;
type "The mean frequency is $(diststats.mean)";
The following example shows how to perform signal filtering using the fft_filters X-Function:
// Import some data with noise and create graph
newbook;
string filepath$ = "Samples\Signal Processing\";
string filename$ = "Signal with High Frequency Noise.dat";
fname$ = system.path.program$ + filepath$ + filename$;
impasc;
plotxy iy:=(1,2) plot:=line;
// Perform low pass filtering
fft_filters filter:=lowpass cutoff:=1.5;
|