2

I am trying to get a real wavelet, w, which is a column vector, from its single-sided Discrete Fourier Transform. According to the theory, the negative frequency side is complex-conjugate of the positive frequency side, but in implementing it in Matlab (using the ifft function) is giving me a headache.

Below, I am listing a small program that transforms a damped sine wavelet, w, into frequency domain as W, then extracting the positive part and augmenting it with conj(flipud(W)), but the inverse FFT of it looks like my input wavelet amplitude modulated with something else . However, w = ifft(W,'symmetric') works fine. Any suggestions to identify the problem will be highly appreciated.

Here is the listing:

clc; clear all
% Genetate a damped sine wavelet
n = 100;
n2 = floor(n/2)+ 1;
dt = .25;
for i  = 1:n
    t = (i-1)*dt;
    w(i,1) = 100 * sin(t) * exp(-0.2*t);
end

figure; subplot(3,2,1);  plot(w);
title('The Signal')
%-------------------------------------
W1  = fft(w);                 % 2-sided
n2 = floor(n/2)+ 1;
W2  = fft(w,n2);              % 1-sided
subplot(3,2,3);plot(real(W2));
title('2-sided abs(W2)')
subplot(3,2,5);plot(imag(W2));
title('2-sided angle(W2)')
%-------------------------------------

w1 = ifft( W1 ) ;                 % Works fine
subplot(3,2,2); plot( w1);
title( ' w2 = ifft(W2);   (2-sided) ' );

% --------------------------------------
% Use the /symmetric/ option of ifft() with
% the single-sided spectrum

w2 = ifft(W2 , 'symmetric');  % 1-sided, works fine

subplot(3,2,4);plot(w2,'k');
title( 'w2 = ifft(W2, "symmetric" )')

% --------------------------------------
% Calculate the complex-cojugate of 1-sided W2
% (excluding the zero frequency point?!), flip it,
% and attach it to the tail of W2 col vector.

H  = flipud(conj(W2(2:n2)));
W3  = [W2 ; H];
w3 = ifft( W3 ) ;    % sourse of my migraine headache
% If you let n =1000 instead of the 100, the effect of
% amplitude-modulation-like effect is less and the output
% (bottom right graph)resembles the input wavelet but
% with a thicker line.
% If n=100 and W2(1:n2-1) in H = ... is used instead
% of the W2(2:n2), you'll get a flying bold eagle!


subplot(3,2,6);plot(w3,'k');
title('w3 = ifft([W2 ; H]')
%---end of the program-------------------
Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
Turhan Tisinli
  • 21
  • 4

3 Answers3

1

The problem is with this line:

W2  = fft(w,n2);              % 1-sided

Unlike the implied assumption that this would return the first n2 outputs of the full length(w)-sized FFT (if that were the case it would give you the expected 1-sided spectrum), this line instead returns the full FFT (2-sided spectrum) of the sequence w truncated to n2 samples.

The fix is then to compute a full FFT of W and then pick the first n2 samples of the result, as your updated code does with:

W1 = fft(w);
W2 = W1(1:n2);
SleuthEye
  • 14,379
  • 2
  • 32
  • 61
0

The following works but I haven't yet figured out why the previous one didn't:

clc; clear all

% Genetate a damped sine wavelet
n = 101;  n2 = floor(n/2) + 1; dt = 0.25;

t = [(0:1:n-1)*dt]'; w = sin(t).*exp(-0.2*t);

figure; subplot(2,1,1);plot(w); 
title('The Wavelet')

W1   = fft(w);                % 2-sided
W2  = W1(1:n2);               % 1-sided

H  = flipud ( conj( W1(2:n2) ) );
W3 = [W2 ; H];  
w3 = ifft(W3);                % 2-sided

subplot(2,1,2); plot(w3); 
title('w3 = ifft( [ W3; H ]')

%---------- end of the program----------
Turhan Tisinli
  • 21
  • 4
0

Down is an augmented version of your sample code at the very begining. The fundamental issue is about the ifft at the end. Consider scaling the real part with (pi./df).

Your 'augmented code' as follows:

close all; clear all; clc; 

% Genetate a damped sine wavelet
n = 512;  n2 = floor(n/2) + 1; dt = 0.25; 
fs=1/dt;   % Digitised data should ever have a sampling rate -right! 

f=linspace(0,fs,n); % Your frequency axis 
df=mean(diff(f)); % Your incrementation on the frequency axis
fc=f(12); % Just to get a monochromatic data in time, frequency is needed    
          % f(12) just to obey the sampling constraints not more than f(n/2+1)


t = [(0:1:n-1)*dt]; w = cos(2.*pi.*fc.*t).*exp(-0.2*t);

figure; subplot(2,1,1);plot(t,w); 
title('The Wavelet')

W1   = fft(w)./n;                % 2-sided
W2  = W1(1:n2);               % 1-sided

H  = fliplr (conj(W1(2:n2-1))); % Intended to stick with the data length...
W3 = [W2 , H];  
w3 = (pi./df).*real(ifft(W3));                % 2-sided (see the scaling)

subplot(2,1,2); plot(w3); 
title('w3 = ifft( [ W3; H ]')

fprintf('Data size at the beginning:\n');
size(w)
fprintf('Final sata size:\n');
size(w3)
Tolga Karabiyikoglu
  • 1