Conversion from FFT to IFFT

Hi,

I have a seismic signal that I have ran thorugh Transfer Functions, with converting it from time-domain to frequency domain using FFT function. Once the transfer functions have been applied, I have tried to bring the signal back to time domain to IFFT, but I am getting only the upper half of the signal (picture attached). Please dont mind the big peak at the end, still need to remove the noise.

G = ifft(Y_corr);

figure

plot(abs(G))

Y_corr is the seismic signal after transfer functions, which needs to be brougth back in time domain. How can I seperate the upper and lower halves of the signal, since "abs" is not really doing that?

Thank You!!!

Answer 1

[Star Strider]

It is likely not efficient to do filtering by editing the fft result (although it can be done). If you want to do it that way, the easiest way is to use the fftshift function to create a symmetrical (positive and negative frequency) version of the Fourier transform. Then, apply the ideal filter (or whatever filter you have available) to both sides of the two-sided Fourier transform. After that, apply ifftshift and then take the ifft of that vector.

The shifts would go something like this —

v = 1:9

v = 1×9

1 2 3 4 5 6 7 8 9

v1 = fftshift(v)

v1 = 1×9

6 7 8 9 1 2 3 4 5

v2 = ifftshift(v1)

v2 = 1×9

1 2 3 4 5 6 7 8 9

Experiment with those functions to understand how they work.

However, that is definitely the long way around. Instead, if you want to pass only a range of frequencies, use the bandpass function. (For best results, include the 'ImpulseResponse','iir' name value pair to design an efficient elliptic filter.) There are other ways to design filters if you have the Signal Processing Toolbox and do not have the bandpass function. I will help with that if necessary.

.

[Nerma Caluk]

What is the exact prupose of fftshift?

This is part of the main code that I currenlty have, which actually gave me some results but I do not know if they are correct:

%% SPZ = data point

%% Fs = sampling frequency

Y = fft(SPZ)

f = Fs*(0:length(Y)/2)/length(Y);

w = 2*pi*f;

%% TfA = Transfer Functions for Acceleration, through which I have ran the omega (w) value

%% Making a irror image of transfer function

Tf1 = [TfA fliplr(TfA)];

Tf1 = Tf1(1:length(Tf1)-1);

Y_corr_A = Y./Tf1'; %% Dividing the FFT values with the transfer function (requirement for this process)

Y_corr_A = Y_corr_A(2:length(Y_corr_A)); %% Eliminating the 1st value of fft since it is -Inf

G_A = ifft(Y_corr_A);

figure

plot(real(G_A))

title('Acceleration vs Time(sec)')

xlabel('Time (sec)')

ylabel('Acceleration (m/s^2)')

The following is the graph plot of abs(Y) and then a plot of real(G_A)

However, if I do apply the fftshift command, as you have instructued I get the following:

V = fftshift(Y);

Y_corr_A = V./Tf1'; % [DU]/[DU/(m/s^2)] = [m/s^2]

Y_corr_A = Y_corr_A(2:length(Y_corr_A));

Y_corr_B = ifftshift(Y_corr_A);

G_A = ifft(Y_corr_B);

figure

plot(real(G_A))

title('Acceleration vs Time(sec)')

xlabel('Time (sec)')

ylabel('Acceleration (m/s^2)')

The ifftshift significanlty increases the values. Does that have to do anything with the mirroring of the Transfer functions that I did?

Thank You so much for your help!

[Star Strider]

‘What is the exact prupose of fftshift?’

It makes the fft symmetric about the origin (D-C or 0 Hz). It is then easier to apply a transfer function to it.

To illustrate this, here is an example of filtering a random signal with a passband of 99 to 100 Hz in the frequency domain, demonstrating the uses of fftshift and ifftshift as well as fft and ifft.

Example —

Fs = 500;
Fn = Fs/2;
L = 5001;
s = randn(1, L);
t = linspace(0, L-1, L)/Fs;

FTs = fft(s)/L;
FTss = fftshift(FTs);
Fv = linspace(-1, 1, L)*Fn;

figure
plot(Fv, abs(FTs))
grid
xlabel('Frequency')
ylabel('Magnitude')

FilterFreq = (Fv <= -99) & (Fv >= -100) | (Fv >= 99) & (Fv <= 100);     % Symmetricalk Passbands

FTss_filtered = FTss;                                                   % Create Filtered Frequency Vector
FTss_filtered(~FilterFreq) = 0;                                         % Set Frequency Components Outside Of The Passband To Zero

figure
plot(Fv, abs(FTss_filtered))
grid
xlabel('Frequency')
ylabel('Magnitude')
title('Filtered Signal')

FTss_filtered = ifftshift(FTss_filtered);
s_filtered = ifft(FTss_filtered);

figure
plot(t, s_filtered)
grid
xlabel('Time')
ylabel('Amplitude')

peridx = islocalmax(s_filtered);
t_peak = t(peridx);
s_filt_freq = 1/mean(diff(t_peak))

s_filt_freq = 99.4791

figure
plot(t, s_filtered)
grid
xlabel('Time')
ylabel('Amplitude')
xlim([1.9  2.1])

The filtered result has a frequency of about 99.5 Hz, corresponding to the desired filter passband frequency of 99 to 100 Hz.

NOTE — I definitely do not recommend this approach! It is much easier to use the time domain filter functions and the filtfilt function (if the filtering function does not also do the filtering) to filter signals, using either FIR or IIR filters.

.

[Walter Roberson]

fft returns a vector of coefficients, first positive frequencies, then negative frequencies. fftshift rearranges to negative then positive. When you apply a frequency filter the assumption of the filter is almost always that the coefficients are in order of increasing frequency, so you fftshift, filter, then undo the fftshift

[Nerma Caluk]

I understand now, thank You both for the clarification!

The reason I cannot apply the transfer functions in time domain is because they are very complex and developed in frequency domain. The process of deconvolution requires that the Magnitudes of the FFT be divided by the given instrument response functions (transfer functions). Due to the fact that I am dividing the amplitudes of FFT with transfer function, I am not sure if fftshift is still applicable. The results give me two different graphs once the fftshift was applied and when not, and both of them are apparenlty larger and off from the assumed, final values, and I do not know what.

Thank You again for the help, if You have any other recommendations, please let me know!

[Star Strider]

That is the best that I can do.

I still do not understand the reason that you cannot do time domain filtering.

[Nerma Caluk]

Thank you for all the help! So if I am to try to solve this problem in time domain, is there a way to convert the frequnecy-domain transfer functions into time-domain? My given transfer function is comprised of multiplication of 4 transfer functions, and I do not know a way to convert them to time-domain in MATLAB.

Thank you again!

[Star Strider]

If you provide the transfer functions, I may be able to convert them into time domain filters (for example using the firls function) that you can use to filter them in the time domain.

It may also be easier than that.

I have not seen them yet, so I do not know if that is possible.

Also, one consiideration with FIR filters is the signal length. How long is (how many elements are in) the signal vector you want to filter?

.

[Nerma Caluk]

The signal I am looking into has 300005 data points, and it is very detailed and heavy. The following are the transfer functions that I have inputed as a function in MATLAB:

function A_SP = transferFunctionA(w)

% Given values for Peak Response

s = w*1j;

R_g = 1800; % coil resistance [ohms]

R_s = 2680; % damping resistance [ohms]

G_1 = 175; % generator constant of the magnet-coil system [V/(m/s^2)]

G_2 = 23700; % pre-amplifier gain

omega_b = 0.31416; % output high pass cut-off angular frequency (rad./s)

omega_p = 57.1199; % output low pass cut-off angular frequency (rad./s)

f_0 = 1; % responant frequency of the pendulum (Hz)

h = 0.85; % damping constant

S_R = 53; % sampling rate, nominal (Hz)

omega_0 = 2*pi*f_0;

K = 204.8;

% Resistance ratio of the damping circuit

G = R_s./(R_g + R_s);

% Transfer function of the seismometer for accleration

S_p = s./(s.^2 + 2.*h.*omega_0.*s + omega_0.^2);

% Transfer function of 8-pole output low-pass anti-aliasing filter

F_l = ([omega_p.^2./(s.^2 + 2.*cos(pi/8).*(omega_p).*s + (omega_p).^2)].^2 ).* [omega_p.^2./(s.^2 + 2.*cos((3*pi)/8).*(omega_p).*s + (omega_p).^2)].^2;

% Transfer function of single-pole high-pass filter in the output amplifier

F_a = s./(omega_b + s);

% Seismometer response for acceleration in flat-response mode

A_SP = K.*G.*G_1.*G_2.*S_p.*F_a.*F_l; %units V./(m./s.^2)

end

Thank You!

[Star Strider]

The ‘w’ argument is obviously the radian frequency. What values should it take?

What is the sampling frequency of the signal? (It should be regular, such that the intervals between samples in the time domain are all equal.)

I decided to let the Control System Toolbox do the ‘heavy lifting’ here —

s = tf('s');            % <— CHANGED
R_g = 1800;             % coil resistance [ohms]
R_s = 2680;             % damping resistance [ohms]
G_1 = 175;              % generator constant of the magnet-coil system [V/(m/s^2)]
G_2 = 23700;            % pre-amplifier gain
omega_b =   0.31416;    % output high pass cut-off angular frequency (rad./s)
omega_p =   57.1199;    % output low pass cut-off angular frequency (rad./s)
f_0 =       1;          % responant frequency of the pendulum (Hz)
h =         0.85;       % damping constant
S_R = 53;               % sampling rate, nominal (Hz)
omega_0 = 2*pi*f_0;
K = 204.8;
% Resistance ratio of the damping circuit 
% G = R_s./(R_g + R_s);
% % Transfer function of the seismometer for accleration
% S_p = s./(s.^2 + 2.*h.*omega_0.*s + omega_0.^2);
% % Transfer function of 8-pole output low-pass anti-aliasing filter
% F_l = ([omega_p.^2./(s.^2 + 2.*cos(pi/8).*(omega_p).*s + (omega_p).^2)].^2 ).* [omega_p.^2./(s.^2 + 2.*cos((3*pi)/8).*(omega_p).*s + (omega_p).^2)].^2;
% % Transfer function of single-pole high-pass filter in the output amplifier
% F_a = s./(omega_b + s);
% % Seismometer response for acceleration in flat-response mode
% A_SP = K.*G.*G_1.*G_2.*S_p.*F_a.*F_l;    %units V./(m./s.^2)

% s = tf('s');
G = R_s/(R_g + R_s);
S_p = s/(s^2 + 2*h*omega_0*s + omega_0^2)

S_p =
 
            s
  ---------------------
  s^2 + 10.68 s + 39.48
 
Continuous-time transfer function.

figure
bodeplot(S_p)
grid

F_l = ([omega_p^2/(s^2 + 2*cos(pi/8)*(omega_p)*s + (omega_p)^2)]^2 )* [omega_p^2/(s^2 + 2*cos((3*pi)/8)*(omega_p)*s + (omega_p)^2)]^2

F_l =
 
                                                       1.133e14
  ------------------------------------------------------------------------------------------------------------------
  s^8 + 298.5 s^7 + 4.456e04 s^6 + 4.299e06 s^5 + 2.908e08 s^4 + 1.403e10 s^3 + 4.743e11 s^2 + 1.037e13 s + 1.133e14
 
Continuous-time transfer function.

figure
bodeplot(F_l)
grid

F_a = s/(omega_b + s)

F_a =
 
      s
  ----------
  s + 0.3142
 
Continuous-time transfer function.

figure
bodeplot(F_a)
grid

A_SP = K*G.*G_1*G_2*S_p*F_a*F_l

A_SP =
 
                                                                            5.758e22 s^2
  -----------------------------------------------------------------------------------------------------------------------------------------------------------------
  s^11 + 309.5 s^10 + 4.788e04 s^9 + 4.802e06 s^8 + 3.399e08 s^7 + 1.741e10 s^6 + 6.411e11 s^5 + 1.619e13 s^4 + 2.478e14 s^3 + 1.696e15 s^2 + 4.982e15 s + 1.405e15
 
Continuous-time transfer function.

figure
bodeplot(A_SP)
grid

Once I get the signal sampling frequency, it will be possible to turn all these into discrete filters.

Meanwhile, you need to determine if these are what you want to use to filter your signal.

.

[Nerma Caluk]

Yes, those all look correct! I have already plotted all of them through 'bode' diagrams, and then used lsim function in one of trys that You have recommended couple of weeks ago, but that function conducted convolution instead of deconvolutuion (dividing) of fft and Trasfer function values.

The sampling frequency is 53Hz. The 'w' values is are radians frequnecy, and should simply be w = 2*pi*f

[Star Strider]

With the sampling frequency provided —

s = tf('s');            % <— CHANGED
R_g = 1800;             % coil resistance [ohms]
R_s = 2680;             % damping resistance [ohms]
G_1 = 175;              % generator constant of the magnet-coil system [V/(m/s^2)]
G_2 = 23700;            % pre-amplifier gain
omega_b =   0.31416;    % output high pass cut-off angular frequency (rad./s)
omega_p =   57.1199;    % output low pass cut-off angular frequency (rad./s)
f_0 =       1;          % responant frequency of the pendulum (Hz)
h =         0.85;       % damping constant
S_R = 53;               % sampling rate, nominal (Hz)
omega_0 = 2*pi*f_0;
K = 204.8;
% Resistance ratio of the damping circuit 
% G = R_s./(R_g + R_s);
% % Transfer function of the seismometer for accleration
% S_p = s./(s.^2 + 2.*h.*omega_0.*s + omega_0.^2);
% % Transfer function of 8-pole output low-pass anti-aliasing filter
% F_l = ([omega_p.^2./(s.^2 + 2.*cos(pi/8).*(omega_p).*s + (omega_p).^2)].^2 ).* [omega_p.^2./(s.^2 + 2.*cos((3*pi)/8).*(omega_p).*s + (omega_p).^2)].^2;
% % Transfer function of single-pole high-pass filter in the output amplifier
% F_a = s./(omega_b + s);
% % Seismometer response for acceleration in flat-response mode
% A_SP = K.*G.*G_1.*G_2.*S_p.*F_a.*F_l;    %units V./(m./s.^2)

% s = tf('s');
Ts = 1/53;                                      % Sampling Period

G = R_s/(R_g + R_s);
S_p = s/(s^2 + 2*h*omega_0*s + omega_0^2);
S_pz = c2d(S_p, Ts)

S_pz =
 
   0.01705 z - 0.01705
  ----------------------
  z^2 - 1.805 z + 0.8175
 
Sample time: 0.018868 seconds
Discrete-time transfer function.

figure
bodeplot(S_p)
grid

F_l = ([omega_p^2/(s^2 + 2*cos(pi/8)*(omega_p)*s + (omega_p)^2)]^2 )* [omega_p^2/(s^2 + 2*cos((3*pi)/8)*(omega_p)*s + (omega_p)^2)]^2;
F_lz = c2d(F_l, Ts)

F_lz =
 
  2.366e-05 z^7 + 0.002952 z^6 + 0.02601 z^5 + 0.04929 z^4 + 0.02647 z^3 + 0.004017 z^2 + 0.0001301 z + 2.964e-07
  ---------------------------------------------------------------------------------------------------------------
        z^8 - 2.794 z^7 + 4.077 z^6 - 3.759 z^5 + 2.352 z^4 - 1.006 z^3 + 0.2832 z^2 - 0.04727 z + 0.00358
 
Sample time: 0.018868 seconds
Discrete-time transfer function.

figure
bodeplot(F_l)
grid

F_a = s/(omega_b + s);
F_az = c2d(F_a, Ts)

F_az =
 
    z - 1
  ----------
  z - 0.9941
 
Sample time: 0.018868 seconds
Discrete-time transfer function.

figure
bodeplot(F_a)
grid

A_SP = K*G.*G_1*G_2*S_p*F_a*F_l;
A_SPz = c2d(A_SP, Ts)

A_SPz =
 
  26.4 z^10 + 7007 z^9 + 9.493e04 z^8 + 1.446e05 z^7 - 2.543e05 z^6 - 2.29e05 z^5 + 1.362e05 z^4 + 9.014e04 z^3 + 1.019e04 z^2 + 212.8 z + 0.2471
  -----------------------------------------------------------------------------------------------------------------------------------------------
      z^11 - 5.593 z^10 + 14.51 z^9 - 23.28 z^8 + 25.79 z^7 - 20.72 z^6 + 12.3 z^5 - 5.38 z^4 + 1.693 z^3 - 0.3636 z^2 + 0.04776 z - 0.002909
 
Sample time: 0.018868 seconds
Discrete-time transfer function.

figure
bodeplot(A_SP)
grid

b = A_SPz.Numerator

b = 1×1 cell array
    {[0 26.3972 7.0068e+03 9.4930e+04 1.4463e+05 -2.5430e+05 -2.2904e+05 1.3620e+05 9.0140e+04 1.0192e+04 212.7577 0.2471]}

a = A_SPz.Denominator

a = 1×1 cell array
    {[1 -5.5932 14.5092 -23.2791 25.7894 -20.7184 12.2972 -5.3797 1.6933 -0.3636 0.0478 -0.0029]}

Use ‘b’ and ‘a’ with the Signal processing Toolbox to filter the signals. It might be appropriate to use the tf2sos funciton to produce a stable filter from them, then use that result with filtfilt.

I do not understand what you intend by deconvolving them. My impression was always that you want to filter the signal with them. In any event, you can use these vectors with the deconv function if necessary, since they are now appropriate for a discrete signal.

.