Changing the inverse fast Fourier transform (ifft) to use an arbitrary waveform instead of sine waves to create a new signal

Question

I know that an inverse fast Fourier transform (ifft) sums multiple sine waves together from data obtain from doing an fft on a signal. Is there a way to create a signal using a new type of inverse fast Fourier transform (ifft) using an arbitrary waveform instead of just using sine waves?

I'm not trying to re-create the original signal. I'm trying to create a new signal using a new type of inverse fast Fourier transform (ifft) using a given arbitrary waveform based on the (frequency, amplitude, phase) data calculated from the fft from the source signal.

The arbitrary waveform is a sampled signal that will replace one period of the sine wave used in the fft. That is, the signal is to be scaled, repeated, and shifted according to the values given by the fft.

See simple example below: the signals I will be applying FFT to are human audio signals about 60 seconds long at 44100 samples (large arrays) so I'm trying to see if I can use / alter the ifft command in some way to create a new signal using / based on an arbitrary waveform.

PS: I'm using Octave 4.0 which is similar to Matlab and the arbitrary waveform signal used to create a new signal will be changed to create different signals.

clear all,clf reset, clc,tic

fs=44100 % Sampling frequency
len_of_sig=2; %length of signal in seconds
t=linspace(0,2*pi*len_of_sig,fs*len_of_sig);
afp=[.5,2.43,pi/9;.3,3,pi/2;.3,4.3,pi/3];  %represents Amplitude,frequency,phase data array

%1 create source signal
ya=0;
for zz=1:size(afp,1)
  ya = ya+afp(zz,1)*sin(afp(zz,2)*t+afp(zz,3));
end

%2 create source frequency domain data
ya_fft = fft(ya);

%3 rebuild original source signal
mag = abs(ya_fft);
phase = unwrap(angle(ya_fft));
ya_newifft=ifft(mag.*exp(i*phase));
ifft_sig_combined_L1=ifft(mag.*exp(i*phase),length(ya_newifft)); 

%4 %%%-----begin create arbitrary waveform to use ---- 
gauss = @(t, t0, g) exp(-((t-t0)/g).^2); % a simple gaussian

t_arbitrary=0:1:44100; % sampling
t_arbitrary_1 = 10000; % pulses peak positions (s)
t_arbitrary_2 = 30000; % pulses peak positions (s)

g = 2000; % pulses width (at 1/e^2) (s)

lilly = gauss(t_arbitrary, t_arbitrary_1, g) - (.57*gauss(t_arbitrary, t_arbitrary_2, g)); %different amplitude peaks
%%%%-----End arbitrary waveform to use---- 

%5 plot
t_sec=t./(2*pi); %converts time in radians to seconds
t_arbitrary_sec=t_arbitrary./length(lilly); %converts time in radians to seconds

subplot(4,1,1);
plot(t_sec,ya,'r')
title('1) source signal')

subplot(4,1,2);
plot(t_sec,ifft_sig_combined_L1)
title('2) rebuilt source signal using ifft')

subplot(4,1,3);
plot(t_arbitrary_sec,lilly,'r')
title('3) arbitrary waveform used to create new signal')

Added a work-flow chart below with simple signals to see if that explains it better:

Section 1) The audio signal is read into an array
Section 2) FFT is done on the signal
Section 3 Red) Normally Inverse FFT uses sin waves to rebuild the signal see signal in red
Section 3 Blue) I want to use an arbitrary signal wave instead to rebuild the signal using the FFT data calculated in (Section 2)
Section 4) New signals created using a new type of Inverse FFT (Section 3).
Please note the new type of Inverse FFT final signal (in blue ) must use the FFT data taken from the original signal.
The signal Sample rate tested should be 44100 and the length of the signal in seconds should be 57.3 seconds long.  I use these numbers to test that the array can handle large amounts and that the code can handle non even numbers in seconds.

Answer 1

Let's start with a function lilly that takes a frequency, an amplitude and a phase (all scalars), as well as a signal length N, and computes a sine wave as expected for the inverse DFT (see note 2 below):

function out = lilly(N,periods,amp,phase)
persistent t
persistent oneperiod
if numel(t)~=N
   disp('recomputung "oneperiod"');
   t = 0:N-1;
   oneperiod = cos(t * 2 * pi / N);
end
p = round(t * periods + phase/(2*pi)*N);
p = mod(p,N) + 1;
out = amp * oneperiod(p);

I have written this function such that it uses a sampled signal representing a single period of the since wave.

The following function uses the lilly function to compute an inverse DFT (see note 1 below):

function out = pseudoifft(ft)
N = length(ft);
half = ceil((N+1)/2);
out = abs(ft(1)) + abs(ft(half)) * ones(1,N);
for k=2:half-1
   out = out + lilly(N,k-1,2*abs(ft(k)),angle(ft(k)));
end
out = out/N;

Now I test to verify that it actually computes the inverse DFT:

>> a=randn(1,256);
>> b=fft(a);
>> c=pseudoifft(b);
recomputung "oneperiod"
>> max(abs(a-c))
ans =  0.059656
>> subplot(2,1,1);plot(a)
>> subplot(2,1,2);plot(c)

The error is relatively large, due to the round function: we're subsampling the signal instead of interpolating. If you need more precision (not likely I think) you should use interp1 instead of indexing using round(p).

Next, we replace the sine in the lilly function with your example signal:

function out = lilly(N,periods,amp,phase)
persistent t
persistent oneperiod
if numel(t)~=N
   disp('recomputung "oneperiod"');
   t = 0:N-1;
   %oneperiod = cos(t * 2 * pi / N);
   gauss = @(t,t0,g) exp(-((t-t0)/g).^2); % a simple gaussian
   t1 = N/4;   % pulses peak positions (s)
   t2 = 3*N/4; % pulses peak positions (s)
   g = N/20;   % pulses width (at 1/e^2) (s)
   oneperiod = gauss(t,t1,g) - (.57*gauss(t,t2,g)); %different amplitude peaks
   oneperiod = circshift(oneperiod,[1,-round(N/4)]); % this will make it look more like cos
end
p = round(t * periods + phase/(2*pi)*N);
p = mod(p,N) + 1;
out = amp * oneperiod(p);

The function pseudoifft now creates a function composed of your basis:

>> c=pseudoifft(b);
recomputung "oneperiod"
>> subplot(2,1,2);plot(c)

Let's look at a simpler input:

>> z=zeros(size(a));
>> z(10)=1;
>> subplot(2,1,1);plot(pseudoifft(z))
>> z(19)=0.2;
>> subplot(2,1,2);plot(pseudoifft(z))

---

Note 1: In your question you specifically ask to use the FFT. The FFT is simply a every efficient way of computing the forward and inverse DFT. The code above computes the inverse DFT in O(n^2), the FFT would compute the same result in O(n log n). Unfortunately, the FFT is an algorithm built on the properties of the complex exponential used in the DFT, and the same algorithm would not be possible if one were to replace that complex exponential with any other function.

Note 2: I use a cosine function in the inverse DFT. It should of course be a complex exponential. But I'm just taking a shortcut assuming that the data being inverse-transformed is conjugate symmetric. This is always the case if the input to the forward transform is real (the output of the inverse transform must be real too, the complex components of two frequencies cancel out because of the conjugate symmetry).

Answer 2

An IFFT is just a way to implement a IDFT. An IDFT is just a weighted sum of sinusoidal waveforms that are integer periodic in aperture.

If you want, you can take almost any DFT or IDFT algorithm or source code, and replace the sin() function with whatever other function you want for waveform synthesis. You can even use different waveforms for different frequencies, or change the synthesis frequencies to non-integer-periodic-in-aperture, if you wish.

Answer 3

The (inverse) Fast Fourier Transform relies on special properties of sinusoidal functions which allow one to move between the time-domain and frequency domain with a much lower computational cost (O(n.log(n))) than for arbitrary waveforms (O(n^2)). If you change the basis waveforms from sinusoids to something else, in general, you can no longer get the computational advantages of the FFT.

In your case, is sounds like you may want to generate a signal that has the same frequency spectrum as your original signal, but not necessarily the same envelope in the time-domain. As I think you've already implemented in your code, the easiest way to do that is simply to change the phase of each frequency-bin in the output of your FFT, and take the inverse FFT. That will generate a time-domain signal that has quite different appearance from your input signal, but will have the same amount of power in each frequency bin.

One subtlety you may need bear in mind is how to change the phases so that the output signal remains real-valued, rather than involving complex numbers. You can think of the Fourier transform as producing a set of positive and negative frequencies. For a real-valued signal, the corresponding positive and negative frequencies must have amplitudes that are complex conjugates of each other. Assuming that your input signal is real-valued, your FFT will already have this property, but you'll need to arrange that the random phases you apply still preserve this relationship between positive and negative frequencies. In practice, this will mean that you only have about half as many random phases to choose - none for zero-frequency bin, and one each for all the positive frequencies (the first n/2 entries) of your FFT, with the other phases being such that the phase at entry k is -1 times that of the entry at (n-k).