matlabfftoctaveifft

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


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')

Plot Updated

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.

Work-flow


Solution

  • 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)
    

    enter image description here

    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)
    

    enter image description here

    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))
    

    enter image description here


    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).