Slow FFT in ksh93 and awk.


 
Thread Tools Search this Thread
Operating Systems OS X (Apple) Slow FFT in ksh93 and awk.
# 1  
Old 08-30-2018
Slow FFT in ksh93 and awk, ARexx and pure awk.

Well i set myself a challenge to have an FFT function using nothing but ksh93 and awk.
It took some serious jiggery pokery and concentration with all the brackets and '$' characters but here is the result. It is RADIX 2 only, but hey, show me another UNIX shell script that does it. It IS SLOW but is ideal for my needs.
It uses ksh's floating point arithmetic and awk to obtain the SINE and COSINE requirements although even awk is not needed for this as typical Taylor series for these could have been created.
This was a learning curve primarily because ksh has a few minor differences to bash; however this is the result.
(For those that know the AMIGA, I am now thinking of porting this to ARexx although it has no SINE, COSINE or SQRT. I have already done SINE and SQRT for ARexx so COSINE is just a simple SINE change.)
This is in the Public Domain, CC0, here and hopefully there will be some interest from outside our community.
The results of the code below:
Code:
Last login: Thu Aug 30 16:06:02 on ttys000
AMIGA:amiga~> cd Desktop/Code/Shell
AMIGA:amiga~/Desktop/Code/Shell> ksh
AMIGA:uw>   
AMIGA:uw> ./FFT_ksh

Real floating point values:
8.5 3.637821186735 -1.70710678119 0.051425210515 1.5 -0.17274555408 -0.292893218805 0.483499156825 0.5 0.483499156845 -0.292893218815 -0.17274555407 1.5 0.05142521051 -1.7071067812 3.637821186735

Imaginary floating point values:
0 -1.93071440555 -0.5 -0.2414680083 -1 0.465638772885 0.5 -1.223607624345 0 1.223607624355 -0.499999999995 -0.46563877289 1 0.24146800831 0.500000000005 1.930714405535

Absolute values of complex numbers:
8.500000 4.118420 1.778820 0.246883 1.802780 0.496649 0.579471 1.315670 0.500000 1.315670 0.579471 0.496649 1.802780 0.246883 1.778820 4.118420 

Done...
AMIGA:uw> exit
AMIGA:amiga~/Desktop/Code/Shell> _

The script:
Code:
#!/bin/ksh
# KSH Version: 'Version AJM 93u+ 2012-08-01'
# From both ${.sh.version} and ${KSH_VERSION}.
#
# Awk Version: 'awk version 20070501'
#
# OSX 10.13.6, default bash terminal running ksh.
# Darwin Barrys-MacBook-Pro.local 17.7.0 Darwin Kernel Version 17.7.0:
# Thu Jun 21 22:53:14 PDT 2018; root:xnu-4570.71.2~1/RELEASE_X86_64 x86_64

# All of this code assumes functions can use global variables, arrays, etc...
# Only constant required, 11 decimal places.
PI=3.14159265359

fft()
{
	# Check for equal length real and imaginary arrays.
	if [ ${#REAL_ARRAY[@]} -ne ${#IMAG_ARRAY[@]} ] || [ ${#REAL_ARRAY[@]} -le 1 ] || [ ${#IMAG_ARRAY[@]} -le 1 ]
	then
		echo "ERROR 1:"
		echo "Number of elements, both real and imaginary, must be the same."
		exit 1
	fi
	# Check for powers of 2.
	N=${#REAL_ARRAY[@]}
	if [ $(( $N&(N-1) )) -ne 0 ]
	then
		echo "ERROR 2:"
		echo "Padding or cropping is required to the nearest power of 2 elements."
		exit 2
	fi
	REAL=()
	IMAG=()
	# For each element into REAL() and IMAG() arrays...
	for ((K=0; K<=$(( $N-1 )); K++))
	do
		SUMREAL=0.0
		SUMIMAG=0.0
		# For each element from REAL_ARRAY() and IMAG_ARRAY() arrays...
		for ((T=0; T<=$(( $N-1 )); T++))
		do
			ANGLE=$(( 2.0*PI*$T.0*$K.0/$N.0 ))
			# Use awk for the COSINE and SINE values...
			COS_ANGLE=$( awk 'BEGIN { printf("%.11f", cos('$ANGLE')); }' )
			SIN_ANGLE=$( awk 'BEGIN { printf("%.11f", sin('$ANGLE')); }' )
			SUMREAL=$(( $SUMREAL+${REAL_ARRAY[$T]}*$COS_ANGLE+${IMAG_ARRAY[$T]}*$SIN_ANGLE ))
			SUMIMAG=$(( $SUMIMAG-${REAL_ARRAY[$T]}*$SIN_ANGLE+${IMAG_ARRAY[$T]}*$COS_ANGLE ))
		done
		REAL[$K]=$SUMREAL
		IMAG[$K]=$SUMIMAG
	done
}

abs()
{
	SQUARE=$(( (${REAL[COUNT]}**2)+(${IMAG[COUNT]}**2) ))
	ABS=$( awk 'BEGIN { printf(sqrt('$SQUARE')); }' )
}

# 16 samples of single cycle square wave, floating point values must always be positive.
#
# REAL_ARRAY=( 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )
# IMAG_ARRAY=( 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )
#
# 8.000000 5.125830 0.000000 1.799950 0.000000 1.202690 0.000000 1.019590
# 0.000000 1.019590 0.000000 1.202690 0.000000 1.799950 0.000000 5.125830

# This has midway '0.5' value padding at the end to bring to the power of 2.
# This is the running default.
#
REAL_ARRAY=( 1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.5 0.5 0.5 )
IMAG_ARRAY=( 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )
#
# 8.500000 4.118420 1.778820 0.246880 1.802780 0.496650 0.579470 1.315670
# 0.500000 1.315670 0.579470 0.496650 1.802780 0.246880 1.778820 4.118420

# Standard test values.
#
# REAL_ARRAY=( 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 )
# IMAG_ARRAY=( 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 )
#
# 4.000000 2.613130 0.000000 1.082390 0.000000 1.082390 0.000000 2.613130

# This assumes all needed variables are global!
fft

# This lot can be discarded...
echo ""
echo "Real floating point values:"
echo "${REAL[@]}"
echo ""
echo "Imaginary floating point values:"
echo "${IMAG[@]}"
echo ""
echo "Absolute values of complex numbers:"

for COUNT in $( seq 0 1 $(( ${#REAL[@]}-1 )) )
do
	abs
	printf "%f " "$ABS"
done
echo ""
echo ""
echo "Done..."

Enjoy all.
Snippets like these are what keep me from going insane.

------ Post updated at 07:31 PM ------

Hi bakunin...

Quote:
Quote:
Originally Posted by wisecracker
I am an AMIGA fanatic still and was searching for a basic FFT routine for Python 1.4.0 for the AMIGA.
My suggestion is to first do it in Assembler. Doing it in Python (or any other interpreted language, except probably FORTH) is like trying to build a race car - out of Lego. And, second, to do it on a DSP. I am not up to date with the latest DSPs any more, but i suggest to take a look at the C66x series from Texas Instruments or the ADSP-SC57x series from Analog Devices. Both are RT-capable floating-point DSPs and from my personal experience (TI32C80 from TI, i can't remember the name of the AD machine) both manufacturers made great procs. You might also try the still available 56k-series by Motorola. I used these (especially the 56001) a lot back in the days.
Thank you for the info but......
I do bizarre things like this to expand my knowledge on the limits of languages, not the limits of HW.
I even got my USB Arduino Diecimila board to talk to my AMIGA with minor HW mods, uploaded to AMINET, using BWBasic purely as a challenge, that's right BWBasic. It has had hundreds of DLs which shows people ARE interested in esoteric stuff like mine.
Real time applications are of no interest at the moment as AudioScope.sh is by definition slow as is its builtin Spec-An using Python 2.7.x with module imports of scipy and scipy.IO. As these ARE dependencies then I have written FFT code that works from Python 2.0.1 to the current 3.7.0. I now have since modified code to work from Python 1.4.0, the only version to work on a stock AMIGA A1200, also up to 3.7.0. Boy how, Python has changed. These mods can be seriously difficult to work around but now I have my building blocks to work with to adapt into AudioScope.sh for AudioScope.sh's Spec-An section without the need for dependencies.
ATM FFT on ksh was a new baby to see how easy it would be, difficult? YES; impossible? NO. I really do get enjoyment out of stretching languages to strange esoteric limits.
Can I do it in AWK? Dunno, I am not good with AWK but ARexx for the AMIGA WILL be serious fun as it does have basic floating point arithmetic but NO SINE, COSINE, SQRT, etc...

Last edited by wisecracker; 09-14-2018 at 07:49 AM.. Reason: (Attached itself to this post.) Now title changed.
This User Gave Thanks to wisecracker For This Post:
# 2  
Old 08-30-2018
I really like FFTs, did quite some more than 30 years ago...

Calling awk (or any external command) twice in a nested loop will inevitably be slow. Why don't you do all of it within one awk script?

Last edited by RudiC; 08-31-2018 at 05:21 AM..
# 3  
Old 08-31-2018
Hi RudiC...

I just love trying near impossible stuff like this. I watch with awe when I see some of you guys code in the shell on things I know little or nothing about so the code above is the way I can get to learn strange subtleties in languages and apply it to things I do know something about.

so two reasons why I haven't tried awk yet.
1) I wanted to see if it was possible to do it in a shell script.
And 2) I don't know enough awk yet and would love to do it, so it IS on the books...

The code is slow but there is no other at all so it might come onto Google's search list bringing people here.
# 4  
Old 09-12-2018
Hi all...

Well guys, it took a while because of the serious limitations of ARexx for a stock AMIGA A1200 but here is a working FFT for it. ARexx has no ARRAY facility but there is a workaround. ;o)
Boy oh boy, ARexx's arithmetic is so convoluted and floating point precision is not the best.
This was a challenge. I could have generated SIN, COS and SQRT in pure ARexx but decided to use an external ARexx math(s) library instead...
It would be just as easy to create SIN, COS and SQRT in ksh too eliminating awk entirely but awk just made it easier...
Just for the record the AMIGA shell/terminal has a large subset of xterm's escape codes.
An exercise in futility? Maybe but serious fun nevertheless...
This is through FS-UAE AMIGA emulation and my real A1200 setup for easy transfer from this Macbook Pro.
ARexx code:
Code:
/* Simple_FFT_DEMO.rexx */

/* This ARexx version needs this dependency: */
/* http://aminet.net/package/util/rexx/RexxMathLib */
CALL ADDLIB('rexxmathlib.library',0,-30,0)

NUMERIC DIGITS 14

SAY ""
SAY "An experimental method to do an FFT in ARexx using the standard test values."
SAY "It requires rexxmathlib.library for the SIN and COS."
SAY "$VER Simple_FFT_DEMO.rexx_(C)06-09-2018_B.Walker_issued_under_GPL2."
SAY ""

/* Create standard test _ARRAYS_ REAL and IMAGINARY. */
/* These MUST be powers of 2 and greater than 2 in size. */
/* Real values, '1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0' */
/* Imag values, '0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0' */
/* Not really ARRAYS but good enough for this DEMO! */
DO N=0 TO 3 BY 1
	REAL_ARRAY.N=1.0
END
DO N=4 TO 7 BY 1
	REAL_ARRAY.N=0.0
END
DO N=0 TO 7 BY 1
	IMAG_ARRAY.N=0.0
END

/* Display them. */
SAY "Input, REAL: IMAG:"
DO N=0 TO 7 BY 1
	SAY "        "||REAL_ARRAY.N||"   "||IMAG_ARRAY.N
END

/* Only constant required, 14 decimal places. */
PI=3.14159265358979 
/* N should be 8 by default in this DEMO. */
/* SAY "Number of elements = "||N||"..." */
N=N-1

CALL FFT

SAY ""
DO N=0 TO 7 BY 1
	SAY "REAL: "||+OUT_REAL.N||",    IMAG: "||+OUT_IMAG.N
END

SAY ""
SAY "Final FFT values to 5 decimal places:"
STR=""
DO N=0 TO 7 BY 1
	CALL ABS_COMPLEX
	NUMERIC DIGITS 14
	IF +ABS <= 0.000000000001
	THEN
		ABS=0
	ENDIF
	NUMERIC DIGITS 6
	STR=STR||+ABS||" "
END
SAY STR
SAY ""
EXIT

/* Create the subroutines required. */
/* FFT subroutine. */
FFT:
DO K=0 TO N BY 1
	SUMREAL=0.0
	SUMIMAG=0.0
	DO T=0 TO N BY 1
		ANGLE=(2.0*PI*T*K)/(N+1)
		SUMREAL=SUMREAL+(REAL_ARRAY.T*COS(ANGLE))+(IMAG_ARRAY.T*SIN(ANGLE))
		SUMIMAG=SUMIMAG-(REAL_ARRAY.T*SIN(ANGLE))+(IMAG_ARRAY.T*COS(ANGLE))
	END
	OUT_REAL.K=SUMREAL
	OUT_IMAG.K=SUMIMAG
END
RETURN
/* FFT subroutine end. */

/* Absolute value of complex number subroutine. */
ABS_COMPLEX:
	ABS=SQRT( (OUT_REAL.N**2)+(OUT_IMAG.N**2) )
RETURN
/* Absolute value of complex number subroutine end. */

# 5  
Old 09-12-2018
People have wanted a good awk fft for years. From this example I'll try and write an awk version. That should help. It has all the facilities needed.
This User Gave Thanks to Corona688 For This Post:
# 6  
Old 09-12-2018
awk version:

Code:
$ cat fft.awk

# Don't feed it  REAL/IMAG/K/T, those are local variables
function fft(REAL_ARRAY, IMAG_ARRAY, K, T, SUMREAL, SUMIMAG,
        ANGLE, COS_ANGLE, SIN_ANGLE)
{
        for(K=1; K in REAL_ARRAY; K++)
        {
                SUMREAL=SUMIMAG=0;

                for(T=1; T in REAL_ARRAY; T++)
                {
                        ANGLE=(2*3.14159*(T-1)*(K-1))/N;
                        COS_ANGLE=cos(ANGLE);
                        SIN_ANGLE=sin(ANGLE);

                        SUMREAL += REAL_ARRAY[T]*COS_ANGLE + IMAG_ARRAY[T]*SIN_ANGLE;
                        SUMIMAG += -REAL_ARRAY[T]*SIN_ANGLE + IMAG_ARRAY[T]*COS_ANGLE;
                }

                REAL[K]=SUMREAL
                IMAG[K]=SUMIMAG
        }
}

BEGIN {
#        REALSTR="1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.5 0.5 0.5"
#        IMAGSTR="0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"

        # Split into arrays.  INDEXES 1-16, not 0-15!
        N=split(REALSTR, REAL_ARRAY, " ");
        split(IMAGSTR, IMAG_ARRAY, " ");
        # Convert strings to integers
        for(X in REAL_ARRAY) REAL_ARRAY[X] += 0.0
        for(X in IMAG_ARRAY) IMAG_ARRAY[X] += 0.0
        # N=16    # Hardcoded 16 elements

        fft(REAL_ARRAY, IMAG_ARRAY);

        printf("Real:\n");
        for(X=1; X in REAL; X++)
                printf(" %.2f", REAL[X]);
        printf("\n");

        printf("Imaginary:\n");
        for(X=1; X in IMAG; X++)
                printf(" %.2f", IMAG[X]);
        printf("\n");

}

$ awk -f fft.awk -v N=16 -v REALSTR="1.0 1.0 1.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.5 0.5 0.5" IMAGSTR="0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0" /dev/null


Real:
 8.50 3.64 -1.71 0.05 1.50 -0.17 -0.29 0.48 0.50 0.48 -0.29 -0.17 1.50 0.05 -1.71 3.64
Imaginary:
 0.00 1.93 0.50 0.24 1.00 -0.47 -0.50 1.22 -0.00 -1.22 0.50 0.47 -1.00 -0.24 -0.50 -1.93

$

It is easily dozens to hundreds of times faster from removing 100% of all external calls.

Last edited by Corona688; 10-04-2018 at 05:24 PM.. Reason: fixing one typo in SUMIMAG equation
This User Gave Thanks to Corona688 For This Post:
# 7  
Old 09-12-2018
Pure awk FFT.

Hi Corona688...

I just knew you would create a translation and get it working on awk. ;o)
You could remove the SIN_ANGLE and COS_ANGLE and put sin(ANGLE) and cos(ANGLE) directly into the two lines below them, eliminating two lines...

Thanks for the translation, I will see if I can eliminate the Python version in AudioScope.sh and attach your name to it...

Another thing is that there is no need to issue the imaginary values for this basic FFT. As long as the real values are actually, cropped or padded to power(s) of 2 then the imaginary values of 0.0 can be internally generated to the same array length as the real values, eliminating the IMAGSTR initialisation on the command line call.

But nice to see an awk version that I can steal... <wink>

Last edited by wisecracker; 09-14-2018 at 07:45 AM.. Reason: Add title...
Login or Register to Ask a Question

Previous Thread | Next Thread

10 More Discussions You Might Find Interesting

1. OS X (Apple)

FFT for Python 2.0.x to 3.7.0.

Hi guys... This is code that was originally designed to work on an upgraded AMIGA A1200 using Python 2.0.x. Unfortunately it broke inside much later versions, NOT because of the print statement/function but other minor subtleties. So this is the final result tested on various machines including... (0 Replies)
Discussion started by: wisecracker
0 Replies

2. AIX

Ksh93/AIX compatibility

Hi everyone ! Im trying to know from wich version of AIX KSH93 is available ? Internet tell me 6.x and 7.x AIX are available, bue what about 5.x ? Is KSH93 available on AIX 5.x ? Is it the same way to manipulate variables as KSH93 on 7.x ? Thanks for your support and have a nice day ! (2 Replies)
Discussion started by: majinfrede
2 Replies

3. OS X (Apple)

FFT for the AMIGA through ksh88 shell.

I don't know if anyone is interested but I have been meddling with FFT for the AMIGA. (Sadly we AMIGAns don't have these luxuries through any scripting language. Below is a Python snippet that uses the builtin 'cmath' module to work with the lowly Python 2.0.1 for the AMIGA. It is part of a... (0 Replies)
Discussion started by: wisecracker
0 Replies

4. UNIX for Advanced & Expert Users

Ksh93 on Linux compatible with ksh93 on AIX

Hi Experts, I have several shell scripts that have been developed on a Linux box for korn ksh93. If we want to run this software on an AIX 6.1 box that runs ksh88 by default can we just change the she-bang line to reference /bin/ksh93 which ships with AIX as its "enhanced shell" to ensure... (6 Replies)
Discussion started by: Keith Turley
6 Replies

5. Shell Programming and Scripting

Making a faster alternative to a slow awk command

Hi, I have a large number of input files with two columns of numbers. For example: 83 1453 99 3255 99 8482 99 7372 83 175 I only wish to retain lines where the numbers fullfil two requirements. E.g: =83 1000<=<=2000 To do this I use the following... (10 Replies)
Discussion started by: s052866
10 Replies

6. Shell Programming and Scripting

cut, sed, awk too slow to retrieve line - other options?

Hi, I have a script that, basically, has two input files of this type: file1 key1=value1_1_1 key2=value1_2_1 key4=value1_4_1 ... file2 key2=value2_2_1 key2=value2_2_2 key3=value2_3_1 key4=value2_4_1 ... My files are 10k lines big each (approx). The keys are strings that don't... (7 Replies)
Discussion started by: fzd
7 Replies

7. Shell Programming and Scripting

The builtin split function in AWK is too slow

I have a text file that contains 4 million lines, each line contains 2 fields(colon as field separator). as shown: 123:444,555,666,777,888,345 233:5444,555,666,777,888,345 623:454,585,664,773,888,345 ...... Here I have to split the second field(can be up to 40,000 fields) by comma into an... (14 Replies)
Discussion started by: kevintse
14 Replies

8. Shell Programming and Scripting

ksh88 or ksh93

Hi all! Does anybody know how can I check if any UNIX installation has implemented ksh88 or ksh93? Thanks in advance. Néstor. (3 Replies)
Discussion started by: Nestor
3 Replies

9. Shell Programming and Scripting

ksh93 deprecation...

Any means of running ksh93 in a ksh88-mode? Might sound odd, but I want/need to restrict U/Win-developed scripts to correspond to the ksh88 version on my Solaris environment(s). Thanks. (2 Replies)
Discussion started by: curleb
2 Replies

10. UNIX for Dummies Questions & Answers

echo is too slow. HELP with Awk

Hello All, Below is a simple script i worte to find the 208th char in a file. If the char = "C" then I re-direct the line to a file called change.txt. If it is not "C" then I re-direct it to a file called delete.txt. My problem is I have a file 0f 500K lines. this script is very slow. I am... (4 Replies)
Discussion started by: eja
4 Replies
Login or Register to Ask a Question