Conversion from FFT to IFFT
Show older comments
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!!!
Answers (1)
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
v1 = fftshift(v)
v2 = ifftshift(v1)
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.
.
11 Comments
Nerma Caluk
on 13 Nov 2022
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!
‘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))
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
on 13 Nov 2022
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
on 14 Nov 2022
Star Strider
on 14 Nov 2022
That is the best that I can do.
I still do not understand the reason that you cannot do time domain filtering.
Nerma Caluk
on 14 Nov 2022
Star Strider
on 14 Nov 2022
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
on 14 Nov 2022
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)
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
figure
bodeplot(F_l)
grid
F_a = s/(omega_b + s)
figure
bodeplot(F_a)
grid
A_SP = K*G.*G_1*G_2*S_p*F_a*F_l
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
on 14 Nov 2022
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)
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)
figure
bodeplot(F_l)
grid
F_a = s/(omega_b + s);
F_az = c2d(F_a, Ts)
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)
figure
bodeplot(A_SP)
grid
b = A_SPz.Numerator
a = A_SPz.Denominator
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.
.
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
on 11 Nov 2022
on 14 Nov 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
