Properties of 2D discrete FFT

Properties of 2D discrete FFT • Fixed sample size • Size of window has to be a base-2 dimension, 32x32, or 64x64 • Periodicity assumption • Particle image is assumed to be periodic • Aliasing • Correlation data is periodic, peak located at area outside of IC window will be found on the opposite side • Bias error • Large shift (less overlapping pattern) leads to smaller and biased peak

Periodicity assumption • Applying FFT to a domain assumes the domain is periodic.

Periodicity assumption • Applying FFT to a domain assumes the domain is periodic. Periodic condition for the domain

Periodicity assumption • Applying FFT to a domain assumes the domain is periodic. Particle image does not satisfy this condition Periodic condition for the domain

Periodicity assumption • Applying FFT to a domain assumes the domain is periodic. Particle image does not satisfy this condition Bias error 1 True peak Measured peak Correlation Periodic condition for the domain Shift (Ds) Bias error

1 1 g f 0 0 0 N 0 N Solution for the bias error 1D example: convolving 1 0 0 -N/2 N/2 Center 0: 1; ±N/2: 1/2 Adjusting correlation coefficients Multiplying weighting factors Reduce bias error Re-scale to [1/2, 1] Image is convolved with itself

correlation Calculating in digital form shift Subpixel interpolation Reason correlation shift correlation • Solution: subpixel interpolation • curve fitting from 3 points near peak • finding the peak in the curve Curve fitting shift

Schemes for subpixel interpolation • Remark: • Gaussian fitting is more commonly used • Parabolic fitting has a divided-by-zero problem when Ri-1 = Ri = Ri+1

Principle Crosscorrelating one IC window with another IC window with known shift Determination of known shift – normal crosscorrelation Implementation Using normal crosscorrelation method to get the velocity map; Detecting bad vectors and replacing with interpolated vectors; Computing crosscorrelation coefficients by shifting another IC window with the displacement determined by the obtained velocity; Advanced techniques- multiple pass interrogation algorithm

IC P-I P-I P-II P-II Advanced techniques- multiple pass interrogation algorithm Normal crosscorrelation multiple pass interrogation

Advanced techniques – hierarchical approach • Principle • Similar to the multiple pass algorithm but the sampling grid system is gradually refined with reducing IC size simultaneously • Advantage • Able to capture complex structure; • Implementation • Run multiple pass algorithms under adaptive grid systems and correspondingly changed IC window

Example of hierarchical approach

Multiple peak detection • Reason • The strongest peak is not always associated with the correct shift, especially in the area of complex structure, e.g., strong vorticity, strong shear • Implementation • Detecting local maximums and keeping them as “peak candidates”; • Compare velocity with surrounding velocities, ruling out the unreasonable “peak candidates” until the searching status is not changed

Example of multiple peak detection Multiple peak detection Single peak detection

Data validation • Reason of the occurrence of bad vectors • Inhomogeneous seeding • Noise • Complex flow structure • Validation method • Direct comparison • Median Filter • Mean Filter

Data validation (cont’d) Direct comparison U6=Ui-1,j+1 U7=Ui,j+1 U8=Ui+1,j+1 U5=Ui-1,j U9=Ui,j U1=Ui+1,j U4=Ui-1,j-1 U3=Ui,j-1 U2=Ui+1,j-1 Median Filter Mean Filter

x3,y3 x4,y4 x0,y0 xc,yc Dy x1,y1 x1,y1 Dx Interpolation • Linear interpolation

Read Tecplot data file • Format of head before data TITLE = "Import0 in Mae513 Velocity vectors [positions in cm] [velocities in cm/s]" VARIABLES = " Position x "," Position y "," Velocity u "," Velocity v " ZONE T="Data", I=34, J=34, F=POINT DT=( SINGLE, SINGLE, SINGLE, SINGLE ) -3.200000000e+001 -3.200000000e+001 0.000000000e+000 0.000000000e+000 0.000000000e+000 -3.200000000e+001 0.000000000e+000 0.000000000e+000 x y u v

Program list (Matlab) FileName = 'Take_01000.dat'; %% TecPlot data file which saves velocity data fid = fopen(FileName,'r'); Line = fgetl( fid ); Line = fgetl( fid ); Line = fgetl( fid ); Dim = sscanf( Line, 'ZONE T="Data", I=%d, J=%d, F=POINT'); Line = fgetl( fid ); nx = Dim(1); %% Grid number in x direction ny = Dim(2); %% Grid number in y direction u = zeros(nx, ny); v = zeros(nx, ny); x = zeros(nx, 1); y = zeros(1, ny); temp = fscanf(fid, '%e'); s = 1; for j = 1:ny for i = 1:nx x(i, 1) = temp(s); y(1, j) = temp(s+1); u(i,j) = temp(s+2); v(i,j) = temp(s+3); s = s + 4; end end fclose(fid); %% now you have x(1:nx, 1), y(1, 1:ny), u(1:nx, 1:ny), v(1:nx,1:ny) %% you can do further processing of x, y, u, v

Program list (Fortran) !! maximum dimension for u and v parameter (mx = 101, my = 101) real x(mx), y(my), u(mx,my), v(mx,my) character*32 FileName /'Take_01000.dat'/ character*80 Line open(100, file=FileName, status = 'old') read(100, *); read(100, *); read(100, '(A80)') Line; read(100, *) !! The line contains the string of !! 'ZONE T="Data", I=%d, J=%d, F=POINT' ie1 = index(Line, 'I=') + 2 ie2 = index(Line, ', J=') - 1 ie3 = index(Line, 'J=') + 2 ie4 = index(Line, ', F=POINT') - 1 read(Line(ie1:ie2), '(I)') nx read(Line(ie3:ie4), '(I)') ny do j=1, ny do i=1, nx read(100, *) x(i), y(j), u(i,j), v(i,j) enddo enddo close(100)

Program list (C/C++) #include <stdio.h> #include <stdlib.h> // maximum dimension for x and y directions #define MX 101 #define MY 101 int main() { char FileName[32]; int nx, ny; char Line[81]; float x[MX], y[MY], u[MY][MX], v[MY][MX]; strcpy(FileName, "Take_01000.dat"); FILE *fp = fopen(FileName, "r"); if ( !fp ) return -1; fgets(Line, 80, fp); fgets(Line, 80, fp); fgets(Line, 80, fp); fgets(Line, 80, fp); // printf("%s", Line); sscanf(Line, "ZONE T=\"Data\", I=%d, J=%d, F=POINT\n", &nx, &ny); printf("nx=%d ny=%d\n", nx, ny); fgets(Line, 80, fp); for (int j=0; j<ny; j++) for (int i=0; i<nx; i++) { fscanf(fp, "%e %e %e %e", &x[i], &y[j], &u[j][i], &v[j][i]); } fclose(fp); }

Properties of 2D discrete FFT