Audio: stft_process: Add conversion to polar format and back

src/audio/mfcc/mfcc_generic.c

@@ -10,6 +10,8 @@
 #include <sof/audio/component.h>
 #include <sof/audio/audio_stream.h>
 #include <sof/math/auditory.h>
+#include <sof/math/icomplex16.h>
+#include <sof/math/icomplex32.h>
 #include <sof/math/matrix.h>
 #include <sof/math/sqrt.h>
 #include <sof/math/trig.h>

src/audio/mfcc/mfcc_hifi4.c

@@ -11,6 +11,8 @@
 #include <sof/audio/component.h>
 #include <sof/audio/audio_stream.h>
 #include <sof/math/auditory.h>
+#include <sof/math/icomplex16.h>
+#include <sof/math/icomplex32.h>
 #include <sof/math/matrix.h>
 #include <sof/math/sqrt.h>
 #include <sof/math/trig.h>

src/audio/mfcc/mfcc_setup.c

@@ -8,6 +8,8 @@
 #include <sof/audio/component.h>
 #include <sof/audio/audio_stream.h>
 #include <sof/math/auditory.h>
+#include <sof/math/icomplex16.h>
+#include <sof/math/icomplex32.h>
 #include <sof/math/trig.h>
 #include <sof/math/window.h>
 #include <sof/trace/trace.h>

src/audio/stft_process/Kconfig

@@ -1,14 +1,30 @@
 # SPDX-License-Identifier: BSD-3-Clause
 
 config COMP_STFT_PROCESS
-	tristate "Template example component"
+	tristate "STFT processing component"
 	default n
 	select MATH_FFT
 	select MATH_32BIT_FFT
 	select MATH_FFT_MULTI
 	help
-	  Select for stft_process component. Reason for existence
-	  is to provide a minimal component example and use as
-	  placeholder in processing pipelines. As example processing
-	  it swaps or reverses the channels when the switch control
-	  is enabled.
+	  Select for stft_process component. STFT acronym means
+	  short term Fourier transform. It converts audio
+	  to multiple FFTs with selected FFT size, hop, and
+	  window function. Possible signal processing can be
+	  done in it in frequency domain for FFTs that is
+	  efficient for more complex signal processing techniques.
+	  The component converts then the frequency domain
+	  version of signal back to normal PCM audio stream
+	  with inverse STFT.
+
+if COMP_STFT_PROCESS
+
+config STFT_PROCESS_MAGNITUDE_PHASE
+	bool "Convert FFTs to polar magnitude and phase"
+	default n
+	help
+	  Select for processing in polar magnitude and phase
+	  domain. Such complex values format is common for
+	  frequency domain signal processing.
+
+endif # COMP_STFT_PROCESS

src/audio/stft_process/stft_process.h

@@ -63,6 +63,7 @@ struct stft_process_buffer {
 struct stft_process_fft {
 	struct icomplex32 *fft_buf; /**< fft_padded_size */
 	struct icomplex32 *fft_out; /**< fft_padded_size */
+	struct ipolar32 *fft_polar;
 	struct fft_multi_plan *fft_plan;
 	struct fft_multi_plan *ifft_plan;
 	int fft_fill_start_idx; /**< Set to 0 for pad left, etc. */

src/audio/stft_process/stft_process_common.c

@@ -5,6 +5,7 @@
 #include <sof/audio/component.h>
 #include <sof/audio/audio_stream.h>
 #include <sof/math/auditory.h>
+#include <sof/math/icomplex32.h>
 #include <sof/math/matrix.h>
 #include <sof/math/sqrt.h>
 #include <sof/math/trig.h>
@@ -97,6 +98,36 @@ static void stft_do_ifft(struct stft_process_state *state, int ch)
 	stft_process_overlap_add_ifft_buffer(state, ch);
 }
 
+#if CONFIG_STFT_PROCESS_MAGNITUDE_PHASE
+static void stft_convert_to_polar(struct stft_process_fft *fft)
+{
+	int i;
+
+	for (i = 0; i < fft->half_fft_size; i++)
+		sofm_icomplex32_to_polar(&fft->fft_out[i], &fft->fft_polar[i]);
+}
+
+static void stft_convert_to_complex(struct stft_process_fft *fft)
+{
+	int i;
+
+	for (i = 0; i < fft->half_fft_size; i++)
+		sofm_ipolar32_to_complex(&fft->fft_polar[i], &fft->fft_out[i]);
+}
+
+static void stft_apply_fft_symmetry(struct stft_process_fft *fft)
+{
+	int i, j, k;
+
+	j = 2 * fft->half_fft_size - 2;
+	for (i = fft->half_fft_size; i < fft->fft_size; i++) {
+		k = j - i;
+		fft->fft_out[i].real = fft->fft_out[k].real;
+		fft->fft_out[i].imag = -fft->fft_out[k].imag;
+	}
+}
+#endif
+
 static void stft_do_fft_ifft(const struct processing_module *mod)
 {
 	struct stft_comp_data *cd = module_get_private_data(mod);
@@ -110,7 +141,12 @@ static void stft_do_fft_ifft(const struct processing_module *mod)
 		if (num_fft) {
 			stft_do_fft(state, ch);
 
-			/* stft_process(state) */
+#if CONFIG_STFT_PROCESS_MAGNITUDE_PHASE
+			/* Convert half-FFT to polar and back, and fix upper part */
+			stft_convert_to_polar(&state->fft);
+			stft_convert_to_complex(&state->fft);
+			stft_apply_fft_symmetry(&state->fft);
+#endif
 
 			stft_do_ifft(state, ch);
 			cd->fft_done = true;

src/audio/stft_process/stft_process_setup.c

@@ -163,6 +163,9 @@ int stft_process_setup(struct processing_module *mod, int max_frames,
 		goto free_fft_buf;
 	}
 
+	/* Share the fft_out buffer for polar format */
+	fft->fft_polar = (struct ipolar32 *)fft->fft_out;
+
 	fft->fft_fill_start_idx = 0; /* From config pad_type */
 
 	/* Setup FFT */

src/audio/tdfb/tdfb_direction.c

@@ -108,7 +108,7 @@ static inline int16_t tdfb_mic_distance_sqrt(int32_t x)
 
 	xs = Q_SHIFT_RND(x, 24, 12);
 	xs = MIN(xs, UINT16_MAX);
-	return sqrt_int16((uint16_t)xs);
+	return sofm_sqrt_int16((uint16_t)xs);
 }
 
 static int16_t max_mic_distance(struct tdfb_comp_data *cd)

src/include/sof/math/fft.h

@@ -11,6 +11,7 @@
 
 #include <sof/audio/module_adapter/module/generic.h>
 #include <sof/audio/format.h>
+#include <sof/math/icomplex32.h>
 #include <sof/common.h>
 #include <stdbool.h>
 #include <stdint.h>
@@ -30,20 +31,6 @@
 #define FFT_SIZE_MAX		1024
 #define FFT_MULTI_COUNT_MAX	3
 
-struct icomplex32 {
-	int32_t real;
-	int32_t imag;
-};
-
-/* Note: the add of packed attribute to icmplex16 would significantly increase
- * processing time of fft_execute_16() so it is not done. The optimized versions of
- * FFT for HiFi will need a different packed data structure vs. generic C.
- */
-struct icomplex16 {
-	int16_t real;
-	int16_t imag;
-};
-
 struct fft_plan {
 	uint32_t size;	/* fft size */
 	uint32_t len;	/* fft length in exponent of 2 */

src/include/sof/math/icomplex16.h

@@ -0,0 +1,84 @@
+// SPDX-License-Identifier: BSD-3-Clause
+//
+// Copyright(c) 2020-2026 Intel Corporation.
+//
+// Author: Amery Song <chao.song@intel.com>
+//	   Keyon Jie <yang.jie@linux.intel.com>
+
+#include <sof/audio/format.h>
+#include <sof/common.h>
+#include <stdint.h>
+
+#ifndef __SOF_ICOMPLEX16_H__
+#define __SOF_ICOMPLEX16_H__
+
+/* Note: the add of packed attribute to icmplex16 would significantly increase
+ * processing time of fft_execute_16() so it is not done. The optimized versions of
+ * FFT for HiFi will need a different packed data structure vs. generic C.
+ *
+ * TODO: Use with care with other than 16-bit FFT internals. Access with intrinsics
+ * will requires packed and aligned data. Currently there is no such usage in SOF.
+ */
+
+/**
+ * struct icomplex16 - Storage for a normal complex number.
+ * @param real The real part in Q1.15 fractional format.
+ * @param imag The imaginary part in Q1.15 fractional format.
+ */
+struct icomplex16 {
+	int16_t real;
+	int16_t imag;
+};
+
+/*
+ * Helpers for 16 bit FFT calculation
+ */
+static inline void icomplex16_add(const struct icomplex16 *in1, const struct icomplex16 *in2,
+				  struct icomplex16 *out)
+{
+	out->real = in1->real + in2->real;
+	out->imag = in1->imag + in2->imag;
+}
+
+static inline void icomplex16_sub(const struct icomplex16 *in1, const struct icomplex16 *in2,
+				  struct icomplex16 *out)
+{
+	out->real = in1->real - in2->real;
+	out->imag = in1->imag - in2->imag;
+}
+
+static inline void icomplex16_mul(const struct icomplex16 *in1, const struct icomplex16 *in2,
+				  struct icomplex16 *out)
+{
+	int32_t real = (int32_t)in1->real * in2->real - (int32_t)in1->imag * in2->imag;
+	int32_t imag = (int32_t)in1->real * in2->imag + (int32_t)in1->imag * in2->real;
+
+	out->real = Q_SHIFT_RND(real, 30, 15);
+	out->imag = Q_SHIFT_RND(imag, 30, 15);
+}
+
+/* complex conjugate */
+static inline void icomplex16_conj(struct icomplex16 *comp)
+{
+	comp->imag = sat_int16(-((int32_t)comp->imag));
+}
+
+/* shift a complex n bits, n > 0: left shift, n < 0: right shift */
+static inline void icomplex16_shift(const struct icomplex16 *input, int16_t n,
+				    struct icomplex16 *output)
+{
+	int n1, n2;
+
+	if (n >= 0) {
+		/* need saturation handling */
+		output->real = sat_int16((int32_t)input->real << n);
+		output->imag = sat_int16((int32_t)input->imag << n);
+	} else {
+		n1 = -n;
+		n2 = 1 << (n1 - 1);
+		output->real = sat_int16(((int32_t)input->real + n2) >> n1);
+		output->imag = sat_int16(((int32_t)input->imag + n2) >> n1);
+	}
+}
+
+#endif /* __SOF_ICOMPLEX16_H__ */

src/math/fft/fft_32.h → src/include/sof/math/icomplex32.h

@@ -1,13 +1,39 @@
 // SPDX-License-Identifier: BSD-3-Clause
 //
-// Copyright(c) 2020-2025 Intel Corporation. All rights reserved.
+// Copyright(c) 2020-2026 Intel Corporation.
 //
 // Author: Amery Song <chao.song@intel.com>
 //	   Keyon Jie <yang.jie@linux.intel.com>
 
 #include <sof/audio/format.h>
-#include <sof/math/fft.h>
+#include <sof/math/exp_fcn.h>
+#include <sof/math/log.h>
+#include <sof/math/trig.h>
 #include <sof/common.h>
+#include <stdint.h>
+
+#ifndef __SOF_ICOMPLEX32_H__
+#define __SOF_ICOMPLEX32_H__
+
+/**
+ * struct icomplex32 - Storage for a normal complex number.
+ * @param real The real part in Q1.31 fractional format.
+ * @param imag The imaginary part in Q1.31 fractional format.
+ */
+struct icomplex32 {
+	int32_t real;
+	int32_t imag;
+};
+
+/**
+ * struct ipolar32 - Storage for complex number in polar format.
+ * @param magnitude The length of vector in Q2.30 format.
+ * @param angle The phase angle of the vector -pi to +pi in Q3.29 format.
+ */
+struct ipolar32 {
+	int32_t magnitude;
+	int32_t angle;
+};
 
 /*
  * These helpers are optimized for FFT calculation only.
@@ -62,3 +88,26 @@ static inline void icomplex32_shift(const struct icomplex32 *input, int32_t n,
 		output->imag = input->imag >> -n;
 	}
 }
+
+/**
+ * sofm_icomplex32_to_polar() - Convert (re, im) complex number to polar.
+ * @param complex Pointer to input complex number in Q1.31 format.
+ * @param polar Pointer to output complex number in Q2.30 format for
+ *		magnitude and Q3.29 for phase angle.
+ *
+ * The function can be used to convert data in-place with same address for
+ * input and output. It can be useful to save scratch memory.
+ */
+void sofm_icomplex32_to_polar(struct icomplex32 *complex, struct ipolar32 *polar);
+
+/**
+ * sofm_ipolar32_to_complex() - Convert complex number from polar to normal (re, im) format.
+ * @param polar Pointer to input complex number in polar format.
+ * @param complex Pointer to output complex number in normal format in Q1.31.
+ *
+ * This function can be used to convert data in-place with same address for input
+ * and output. It can be useful to save scratch memory.
+ */
+void sofm_ipolar32_to_complex(struct ipolar32 *polar, struct icomplex32 *complex);
+
+#endif /* __SOF_ICOMPLEX32_H__ */

src/include/sof/math/sqrt.h

@@ -1,15 +1,32 @@
 /* SPDX-License-Identifier: BSD-3-Clause
  *
- * Copyright(c) 2021 Intel Corporation. All rights reserved.
+ * Copyright(c) 2021-2026 Intel Corporation.
  *
  * Author: Shriram Shastry <malladi.sastry@linux.intel.com>
  *
  */
 
-#ifndef __SOF_MATH__SQRTLOOKUP__H
-#define __SOF_MATH__SQRTLOOKUP__H
+#ifndef __SOF_MATH_SQRT_H__
+#define __SOF_MATH_SQRT_H__
 
 #include <stdint.h>
 
-uint16_t sqrt_int16(uint16_t u);
-#endif
+/**
+ * sofm_sqrt_int16() - Calculate 16-bit fractional square root function.
+ * @param u Input value in Q4.12 format, from 0 to 16.0.
+ * @return Calculated square root of n in Q4.12 format, from 0 to 4.0.
+ */
+uint16_t sofm_sqrt_int16(uint16_t u);
+
+/**
+ * sofm_sqrt_int32() - Calculate 32-bit fractional square root function.
+ * @param n Input value in Q2.30 format, from 0 to 2.0.
+ * @return Calculated square root of n in Q2.30 format.
+ *
+ * The input range of square root function is matched with Q1.31
+ * complex numbers range where the magnitude squared can be to 2.0.
+ * The function returns zero for non-positive input values.
+ */
+int32_t sofm_sqrt_int32(int32_t n);
+
+#endif /* __SOF_MATH_SQRT_H__ */

src/math/CMakeLists.txt

@@ -16,7 +16,7 @@ if(CONFIG_MATH_LUT_SINE_FIXED)
 endif()
 
 if(CONFIG_SQRT_FIXED)
-  list(APPEND base_files sqrt_int16.c)
+  list(APPEND base_files sqrt_int16.c sqrt_int32.c)
 endif()
 
 if(CONFIG_MATH_EXP)
@@ -91,6 +91,10 @@ if(CONFIG_MATH_MU_LAW_CODEC)
   list(APPEND base_files mu_law.c)
 endif()
 
+if(CONFIG_MATH_COMPLEX)
+  list(APPEND base_files complex.c)
+endif()
+
 is_zephyr(zephyr)
 if(zephyr) ###  Zephyr ###