MATLAB: Time to frequency domain and back again

fftfrequency domainiffttime domain

I am trying to find the proper way to convert a periodic signal from time-domain to frequency-domain (with the proper amplitude and phase) and then retrieving the data again into time-domain. The purpose of this is to do some post processing in frequency-domain and observe the modified signal in time-domain. I am using fft() and ifft() to do so. I tested this approach by creating an arbitrary periodic function consisting of the sum two arbitrary sinusoidal waves (arbitrary phase, amplitude and frequency) and a random noise signal. However, I always get this nonlinear shift in the frequency of the retrieved signal, that is bigger at higher frequency values. The arbitrary signal:

Original signal- retrieved signal comparison (Time-domain)

Original signal- retrieved signal comparison (Frequency domain)

I am trying to investigate the source of this error. I made sure that I am compensating for the windowing effect of the fft as shown in the myfft and myifft functions. I also investigated the effect of frequency step (Spectral leakage) by sampling the total length of the time signal but it didn't seem to affect the results.

The following is the testing code and the used functions:

clc                                                  %Clear display
clear all                                            %Clear all variable
close all                                            %Close all windows
f1=1000;                                             %Frequency of the first sine wave
f2=5000;                                             %Frequency of the second sine wave
Tau=1/max([f1,f2])/100;                              %Sampling time step to make Fsampling = 10*Fmax
Lsignal=[10 100 1000 10^4];                     %Multiplication factor applied to define total signal length to compare the effect of freqnecy resolution
Lsignal=2.^nextpow2(Lsignal);                     %multiplication factor applied to define total signal length
Linecolor=['b','m','k','b','r'];                  %Line colors array used for each total signal length/ freqnecy resolution
for NLsignal=1:length(Lsignal)
    %Create time domain signal
    Ttotal=max([1/f1,1/f2])*Lsignal(NLsignal);       %Total sample time =Lsignal*Longest periodic time of the signal. This value was selected to insure minmum lekage effect, without windowing
    t=0:Tau:Ttotal;                                  %Time vector size is extended
    A1=2; A2=5;                                      %Signals amplitude
    Phi1=135; Phi2=90;                               %Signals Phase shifts in degrees
      y1=2*sin(2*pi*f1*t-Phi1*pi/180);                 %Arbitrary sine wave 1
      y2=5*sin(2*pi*f2*t-Phi2*pi/180);                 %Arbitrary sine wave 2
      y3=0.7*(rand(length(y2),1)'-0.5);                %Arbitary random signal
      yt=y1+y2+y3;                                     %Total sine wave
      %Display orginal time domain signal
      if NLsignal==1
          figure (1)
          plot(t,y1,'color','r','linewidth',2)
          hold on
          plot(t,y2,'k','linewidth',2)
          hold on
          plot(t,y3,'m','linewidth',2)
          hold on
          plot(t,yt,'b','linewidth',2)
          xlabel('Time [sec]')
          ylabel('Amplitude')
          legend ({'Arbitrary sine wave 1';'Arbitrary sine wave 2';'Noise' ;'Total sine wave'});
          xlim([0 max([1/f1,1/f2])])
          ylim([-8 8])
          Fsampling=1/t(2);
          WT=6;                                     % Window type  refer to myFFT help for more info

          title('Single period of the time doamin of the singal')
      else
      end
      %Display orginal time domain for comparison
      if NLsignal==1
          figure (2)
          plot(t,yt,'color',Linecolor(NLsignal),'linewidth',2)
          hold on
          xlabel('Time [sec]')
          ylabel('Amplitude')
          legend ({'Total sine wave'});
          xlim([0 max([1/f1,1/f2])])
          ylim([-8 8])
          Fsampling=1/t(2);
          WT=6;                                     % Window type  refer to myFFT help for more info
          title('Single period of the time doamin of the singal')
      else
      end
      %To frequency domain
      [Yxx,ff] = myfft(Fsampling,yt'-mean(yt),WT);
      figure(3)
      plot(ff,abs(Yxx),'color',Linecolor(NLsignal));
      hold on
      title('Frequency domain signal')
      xlabel('Frequency [Hz]')
      ylabel('Amplitude')
      figure(4)
      plot(ff,angle(Yxx)*180/pi,'color',Linecolor(NLsignal),'linewidth',2)
      xlabel('Frequency [Hz]')
      ylabel('angle [deg]')
      %To time domain
      [ yi ] = myifft((Yxx),length(t),WT);
      %Display Retrieved time domain for comparison
      figure(2)
      hold on
      plot(t,yi,'color',Linecolor(NLsignal),'linestyle','--','linewidth',2);
      %To frequency domain again for verification
      [Yxx2,ff2] = myfft(Fsampling,yi,WT);
      %Display Retrieved frequency domain for comparison
      figure(3)
      hold on
      plot(ff2,abs(Yxx2),'color',Linecolor(NLsignal),'linestyle','--','linewidth',2);
      xlim([0 max([f1,f2])*1.5])
      ylim([0 8])
      figure(4)
      hold on
      plot(ff,angle(Yxx)*180/pi,'color',Linecolor(NLsignal),'linestyle','--','linewidth',2);
      xlim([min([f1,f2])-1 min([f1,f2])+1])
      xlabel('Frequency [Hz]')
      ylabel('angle [deg]')
end
figure(2)
grid on
legend({'Original signal';'Retrieved signal:10*TimePeriod ';'Retrieved signal:10^2*TimePeriod ';'Retrieved signal:10^3*TimePeriod ';'Retrieved signal:10^4*TimePeriod ';'Retrieved signal:10^5*TimePeriod '})
figure(3)
grid on
legend({'10*TimePeriod';'Retrieved signal:10*TimePeriod ';'10^2*TimePeriod';'Retrieved signal:10^2*TimePeriod ';'10^3*TimePeriod';'Retrieved signal:10^3*TimePeriod ';'10^4*TimePeriod';'Retrieved signal:10^4*TimePeriod ';'10^5*TimePeriod';'Retrieved signal:10^5*TimePeriod '})
figure(4)
grid on
legend({'10*TimePeriod';'Retrieved signal:10*TimePeriod ';'10^2*TimePeriod';'Retrieved signal:10^2*TimePeriod ';'10^3*TimePeriod';'Retrieved signal:10^3*TimePeriod ';'10^4*TimePeriod';'Retrieved signal:10^4*TimePeriod ';'10^5*TimePeriod';'Retrieved signal:10^5*TimePeriod '})

myfft function:

function [ fftres, freq ] = myfft( Fs,y,WT)
% [ fftres, freq ] = myfft( Fs,y,WT)

% Inputs

% ---------------------------------------------------------------------------------------------------



% Fs: The original sampling frequency of the acquired signal
% y : Input signal in time domain
% WT : Window type. It is used to define the type of window and apply the corresponding CG correction

% 1=Rectangular window

% 2=Hamming window

% 3=Hanning window

% 4=Bartlett window

% 5=Blackman-Harris window

% 6=Flat Top window

%

% Outputs

% ---------------------------------------------------------------------------------------------------
% fftres: the resulted complex fft
% freq: Corresponding frequency vector
Windowtype = {'rectwin'; 'hamming' ;'hanning';'bartlett';'blackmanharris';'flattopwin'};
%Windowing correction  factors [1] F. J. Harris, “On the use of windows for harmonic analysis with the discrete fourier transform,”Proceedings of the IEEE, vol. 66, no. 1, pp. 51–83, Jan. 1978.

%CG=1.0000; NG=1.0000; %Rectangular window

%CG=0.5400; NG=0.3974;  %Hamming window

%CG=0.5000; NG=0.3750; %Hanning window

%CG=0.5000; NG=0.3333; %Bartlett window

%CG=0.3587; NG=0.2580; %Blackman-Harris window

%CG=0.2156; NG=0.1752; %Flat Top window

CG=[1.0000 0.5400 0.5000 0.5000 0.3587 0.2156];
NG=[1.0000 0.3974 0.3750 0.3333 0.2580 0.1752];
L = length(y);
NFFT = 2^nextpow2(L); 
eval(['Y = fft(y.*',char(Windowtype(WT)),'(L),NFFT)/L;'])
f = Fs/2*linspace(0,1,NFFT/2+1);
fftres = 2*(Y(1:NFFT/2+1))/CG(WT);
freq   = f(1:NFFT/2+1) ;
end
*myifft* function:
function [ Pt ] = myifft( Pf,Lt,WT)
% [ fftres, freq ] = myfft( Fs,y,WT)
% Inputs
% ---------------------------------------------------------------------------------------------------
% Pf: Vector representing frequency-domain complex amplitudes
% Lt : Length of the vector
% WT : Window type. It is used to define the type of window and apply the corresponding CG correction
% 1=Rectangular window
% 2=Hamming window
% 3=Hanning window
% 4=Bartlett window
% 5=Blackman-Harris window
% 6=Flat Top window
%
% Outputs
% ---------------------------------------------------------------------------------------------------
% Pt: the resulted time domain signal
Windowtype = {'rectwin'; 'hamming' ;'hanning';'bartlett';'blackmanharris';'flattopwin'};
%Windowing correction  factors [1] F. J. Harris, “On the use of windows for harmonic analysis with the discrete fourier transform,”Proceedings of the IEEE, vol. 66, no. 1, pp. 51–83, Jan. 1978.
%CG=1.0000; NG=1.0000; %Rectangular window
%CG=0.5400; NG=0.3974;  %Hamming window
%CG=0.5000; NG=0.3750; %Hanning window
%CG=0.5000; NG=0.3333; %Bartlett window
%CG=0.3587; NG=0.2580; %Blackman-Harris window
%CG=0.2156; NG=0.1752; %Flat Top window
CG=[1.0000 0.5400 0.5000 0.5000 0.3587 0.2156];
NG=[1.0000 0.3974 0.3750 0.3333 0.2580 0.1752];
Nfft=length(Pf);
Pt = ifft(Pf,Lt,'symmetric');
LPt=length(Pt);
eval(['Winl=' ,char(Windowtype(WT)),'(LPt);']);
Pt=Pt*LPt/2./Winl*CG(WT);
end

Best Answer

Hi weam,

This kind of error occurs all the time. Because of the frequency shift I was interested to see if you had used nextpow(2) anywhere in your calculations. You had, and after replacing

NFFT = 2^nextpow2(L);    with     NFFT = L;

it works like a charm.

Nextpow can be made to work, although zerofilling has disadvantages in any case. You had the right idea to use NFFT instead of L to construct the frequency grid, so it is not very clear how the error is propagating. But rather to chase that down, it much better and easier besides, to forget about nextpow altogether. You have 10^4 points, so the fft is super fast anyway and there is no need for the complications brought in by going to 2^n points.

At one time 2^n was integral to a breakthrough in computing, but since then ffts have improved and there is no need to use 2^n all the time. It made sense 40 years ago, but it's really time to move on. Unless you have a very large number of points and it is a prime number, there is no need to do any zerofilling.

Related Solutions

MATLAB: How to FFT a portion of an audio signal using window size 8000.

hello

see demo code below, adapt it to your own needs

clc
clearvars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
[signal,Fs] = audioread('test_voice.wav');
[samples,channels] = size(signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 256;    % 
OVERLAP = 0.75;
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 4;
if decim>1
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    Fs = Fs/decim;
    end
   signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;
%%%%%% legend structure %%%%%%%%
for ck = 1:channels
    leg_str{ck} = ['Channel ' num2str(ck) ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),plot(time,signal);grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend(leg_str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);
% apply A weigthing if needed

if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end
figure(2),plot(freq,sensor_spectrum_dB);grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
legend(leg_str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
    % FFT normalisation and conversion amplitude from linear to dB (peak)
    sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
    % apply A weigthing if needed
    if option_w == 1
        pondA_dB = pondA_function(fsg);
        sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
        my_title = ('Spectrogram (dB (A))');
    else
        my_title = ('Spectrogram (dB (L))');
    end
    % saturation of the dB range : 
    % saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
    min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
    sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
    % plots spectrogram
    figure(2+ck);
    imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');
end
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows
    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

Best Answer

Related Solutions

MATLAB: How to FFT a portion of an audio signal using window size 8000.

Related Question