CONVERT FREQUENCY DOMAIN DATA TO TIME DOMAIIN

i have .csv file which contains frequency (Hz), magnitude(dB), phase(degrees). these data are range from 20hz to 2000000Hz(2MHz) and these are non uniformly spaced. i want to convert these data in time domain.what should i do? i have attached the frequency data.

% === STEP 1: Load SFRA Data ===
data = readmatrix('frequency data.csv');       % FAR END SHORTED
freq = data(:,1);                         % Frequency in Hz
mag_dB = data(:,2);                       % Magnitude in dB
phase_deg = data(:,3);                    % Phase in degrees

% === STEP 2: Convert to Linear Magnitude and Radians ===
mag_linear = 10.^(mag_dB / 20);           % Convert dB to linear scale
phase_rad = deg2rad(phase_deg);           % Convert degrees to radians

% === STEP 3: Interpolate to Uniform Frequency Grid ===
N = 2048;                                 % Power-of-2 for IFFT
f_uniform = linspace(min(freq), max(freq), N);
mag_uniform = interp1(freq, mag_linear, f_uniform, 'linear');
phase_uniform = interp1(freq, phase_rad, f_uniform, 'linear');

% === STEP 4: Construct Complex Frequency Response ===
H_f = mag_uniform .* exp(1j * phase_uniform);

% === STEP 5: Create Hermitian Symmetry for Real IFFT ===
H_full = [H_f, fliplr(conj(H_f(2:end-1)))];

% === STEP 6: Apply IFFT ===
h_t = ifft(H_full, 'symmetric');          % Time-domain impulse response

% === STEP 7: Time Vector ===
Fs = 2 * max(f_uniform);                  % Sampling frequency (Nyquist)
t = (0:length(H_full)-1) / Fs;            % Time in seconds

[Star Strider]

You can calculate the inverse Fourier transform of your frequency domain data, however I doubt that the result would be at all meaningful. There is simply too much missing information in your original data.

Answer 1

[William Rose]

@ridham,

I think your approach will work. You convert dB to magn, you convert degrees to radians. You convert magn, phase to real, imag. This is all good. You interpolate the unevenly spaced freqs to a linear set of freqs. You make sure the upper half of the spectrum is the conjugate mirror image of the lower half. You check if the number of freqs is even or odd and adjust as necessary. This is good. The factor of 8 which you use when allocating the vector for the complex vector must depend on your knowledge of the unevenly spaced frequencies. Since we don't have your actual data, we cannot be sure it will all work, but it looks good to me. Test it with simulated data for which you know the correct answer.

Good luck with your work.

[Paul]

Hi William,

To the extent this approach might work at all (hard to say w/o seeing the actualy data as you said), my suggestions would be:

interpolate real and imaginary parts instead of mag and phase. If interpolating in mag and phase, at least make sure the phase is correctly unwrapped before interpolation.

I'd make sure to fill in Hfft for frequency points from 0 to fmin in H_fft_pos before going to Hfull.

Definitely test the approach with simulated data, hopefully simulated data that has characteristics similar to the actual data.

Answer 2

[Paul]

Here's an example that illustrates the potential utility and pitfalls of this approach.

EDIT: 19 Jan 2026, updated code to ensure the extrapolated transform is real at f = 0

Generate time domain data

N = 100;
n = 0:N-1;
xval = (0.94.^n).*(n>=10 & n <= 30);

Plot its transform (magnitude) so we get an idea of what we are dealing with; note the large magnitude at low frequencies

Xval = fft(xval);
figure
plot(n/N,abs(Xval))
xlabel('freq');ylabel('abs');

Define the range of frequencies over which we'll have non-uniformly spaced samples of frequency domain data. For simplicity, we assume the max frequency is the Nyquist frequency.

% Nyquist frequency
fNyq = 0.5;
fmin = 0;
fmax = fNyq;

Generate frequency domain samples over random frequencies between fmin and fmax. We use N frequencies, which is what would be needed to represent the time domain signal if they were uniformly spaced.

rng(101);
fsamp = sort([0,fmax*rand(1,N-2),fmax]);

Evaluate the transform at the frequency domain samples.

Xsamp = freqz(xval,1,fsamp,1);

Overlay the random samples on the transform

hold on
plot(fsamp,abs(Xsamp),'o'),grid

Do the interpolation/extraploation and then transform back to time domain. Here, nmin defines the number of low frequency samples that we exclude (the OP stated that the data is not available down to DC). If nmin = 0, we're just interpolating fsamp as given. If nmin = 1, then we're interpoationg from fsamp(2:end) and extrapolating back to f = 0, and so on.

for nmin = 0:3
helper(n,xval,fsamp(nmin+1:end),Xsamp(nmin+1:end),nmin,fmax);
end

●

We see that for this data set interpolating in the frequency domain over 0-Fnyq works o.k., sort of. But as we start neglecting (high amplitude) lower frequency samples and have to extapolate down to DC the result starts to break down.

I guess the approach proposed by the OP (as modified below) might work depending on the underlying physics, but I can definitely see it being problematic. And I'm sure this all depends on how well the original frequency-domain samples capture the important parts of the transform, among who knows what other things.

Should also be able to use nufft w/o interpolation/extrapolation instead of ifft w/ interpolation/extrapolation for this type of problem.

function helper(n,xval,fsamp,Xsamp,nmin,fmax)

% define number of uniformly spaced frequency domain samples from 0 to
% sampling frequency
N2 = 512;
% define "lower half" frequency vector with uniform spacing with N2 being even
finterp = linspace(0,fmax,N2/2+1);
% linear interpolation/extrapolation of real and imaginary parts
Xinterp = interp1(fsamp,Xsamp,finterp,'linear','extrap');

% Ensure points at Nyquist and DC are real, there may be better ways to do
% this depending on data?
Xinterp(end) = real(Xinterp(end));
Xinterp(1) = abs(Xinterp(1));

% inverse transform
xinterp = ifft(Xinterp,N2,'symmetric');
% same as
% xinterp = ifft([Xinterp,fliplr(conj(Xinterp(2:end-1)))],'symmetric');

% plot
figure
n2 = 0:N2-1;
plot(n,xval,n2,xinterp),grid
legend('original','after interp/extrapolation');
title("nmin = " + double(nmin));
xlabel('n')

end

[ridham]

this is the frequency data i have. please kindly guide me for the code to get correct time domain response.

[Paul]

Let's take a look at the frequency-domain data

data = readmatrix('frequency data.csv');
    
% Extract columns (adjust indices if your CSV has different column order)
freq = data(:, 1);        % Frequency in Hz
mag_dB = data(:, 2);      % Magnitude in dB
phase_deg = data(:, 3);   % Phase in degrees

figure
plot(subplot(211),freq,mag_dB),grid
plot(subplot(212),freq,phase_deg),grid

figure
plot(diff(freq)),grid

It looks very smooth and dense, so interpolation might not be so bad. It looks flat at low frequency with zero phase, so instead of extrapolating to 0, maybe it would be better to prepend a point to the data at f = 0 with value equal to the magnitude of the first data point before interpolation.

Not sure what to make of the phase being so far away from zero (or 180) at the last data point. Maybe someone else has an idea.

Do you have any details on how the frequency-domain data was generated?

[Paul]

Here's an example with a known transfer function that has similar characteristics to get an idea of how one might attack this problem and where some problems might lie.

Define the transfer function

w0 = 50;
H = tf([1,2*.1*w0,w0^2],[1,2*.9*w0,w0^2])*tf(1,[1/1000,1]);

Exact impulse response and a function to evaluate it

ximpulse = ilaplace(poly2sym(H.num{1})/poly2sym(H.den{1}));
ximpfun  = matlabFunction(ximpulse);

Get the continuous-time frequency response and plot it

[m,p,w] = bode(H);
m = squeeze(m); p = squeeze(p);
figure
semilogx(subplot(211),w,db(m),'-o'),grid
semilogx(subplot(212),w,p,'-o'),grid

w = 1e5 rad/sec is high enough to get to "steady state" in the frequency domain where the phase is constant.

bode does not return linearly equally spaced points in frequency, so we interpolate the real and imaginary parts, though linearly interpolating db(gain) and phase in log10 frequency space should work as well.

Compute the complex response

h = m.*exp(1j*p*pi/180);

Number of points to sample in the frequency domain

N = 2^16;

Equally spaced frequency vector and corresponding time samples

wsamp = linspace(0,w(end),N/2+1); % lower half
dw = wsamp(2);
df = dw/2/pi;
Fs = df*N;
Ts = 1/Fs;
t = (0:N-1)*Ts;

Interpolate real and imaginary, ensure conjugate symmetry.

hsamp = interp1([0;w],[m(1);h],wsamp);
hsamp(end) = real(hsamp(end)); % ensure symmetry

Inverse transform back to time domain

xrecon = ifft(hsamp,N,'symmetric')/Ts;

Compare

figure
fplot(subplot(311),ximpulse,[0,1]),grid
hold on
plot(subplot(311),t,xrecon)
xlim('padded')
legend('x','xrecon');
copyobj(get(gca,'children'),subplot(312)),grid
xlim([0,.2]);
copyobj(get(gca,'children'),subplot(313)),grid
xlim([-0.001,0.005]);

The top plot looks pretty good, but there is a small transient at the end. The bottom plot shows that the initial value is 1/2 of what it should be, and that xrecon includes very high frequency oscillation (which is very evident after zooming in further).

I think that these anomalies arise, at least in part, because we are taking the IDFT of frequency-domain samples of a continuous-time transform instead of samples of the corresponding discrete-time transform.

The samples of the true, discrete-time transform can be obtained from time-domain samples of the impulse response

X = fft(ximpfun(t));

Compare in the frequency domain

ii = wsamp < pi/Ts;
figure
semilogx(subplot(211),w,db(m),wsamp(ii),db(abs(Ts*X(ii))),'x'),grid
semilogx(subplot(212),w,p,wsamp(ii),180/pi*angle(X(ii)),'x'),grid
xlabel('rad/sec')

Obviously, if we take the IDFT of the red samples, we'd recover time domain samples of x.

The frequency response of the OP appears to not have reached a "steady state" in the frequency domain (phase a multiple of 90 if continous-time, phase either 0 or 180 if discrete-time). We can pretend it is, and proceed as above. In this example, one could explore what happens if, for example, w(end) = 1e3 rather than 1e5 and then proceed as shown above (I tried that, it didn't look too good, but it might have been good enough). Or one could set w(end) = 1e3 and then try to come up with an extrapolation scheme in the frequency domain to "steady state" perhaps based on known physical characteristics of the system under test. After that extrapolation, one could also try to compute the discrete-time transform via periodic summation of the continuous-time transform. Also, it's not clear that the response of the OP extends to low enough frequency either, though it looks close.

Might want to ask this question on other forums that are more focused on DSP.

Answer 3

[Arkadiy Turevskiy]

Great question and discussion. I wanted to add my take. You can do it all using System Identification Toolbox. It can estimate an LTI model from frequency domain data.

See attached script.

Here is how I did it.

1 Download your csv file. Drag and drop it into MATLAB to open Import tool and import the data.

2. Click "Import Selection" to import. Note that in the drop down menu under import, you can also generate MATLAB code to add to your script so you can do it programmatically.

3. I now have the data as a variable in my workspace. I want to convert it to an frd object that I can use in System Identification Toolbox. I asked MATLAB Copilot to help with that. It provides the code that I can copy to my script and run - awesome!

4. I now have an frd object and can do a bode plot.

5. Now I am ready for system identification. Launch System Identification app from app gallery or type >>systemIdentification. Load sys as Frequency domain, frd object data.

6. Use the app to identify an LTI model that we can use for time-domain analysis. I chose to estimate a state-space model because for state-space the tool automatically suggests optimal model order.

7. The bode plot comparison looks good.

8. Again, as with Import tool, I can get MATLAB code used by the tool so I can do everything programmatically.

9. I can now bring my estimated LTI model to MATLAB workspace and use it for time-domain analysis - bring it to Simulink, design a controller, or just do a step response plot.

I packaged everything I did into the attached script.

Summary:

• Use System Identification Toolbox to convert frequency response data to an LTI model to use for time-domain analysis,
• You can do most steps interactively and as you do that you can also save the code for programmatic reuse.
• Use MATLAB Copilot!

Hope this helps.

Arkadiy

[Paul]

Hi Arkadiy.

I tried running the script but the Bode plot isn't working. The original call to bode in the script didn't work and the call didn't look right to me, so I changed it to bodeplot with what I thought would be the correct arguments but it still doesn't work. Can you correct the code below? I'd like to see what that plot looks like.

%% Set up the Import Options and import the data
opts = delimitedTextImportOptions("NumVariables", 3, "Encoding", "UTF-8");

% Specify range and delimiter
opts.DataLines = [2, Inf];
opts.Delimiter = ",";

% Specify column names and types
opts.VariableNames = ["Frequency", "x_Magnitude_dB__", "x_Degree____"];
opts.VariableTypes = ["double", "double", "double"];

% Specify file level properties
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";

% Import the data
frequencyData = readtable("frequency data.csv", opts);


%% Clear temporary variables
clear opts% Extract columns (by name or by index)
f_Hz      = frequencyData{:,1};   % Hz
mag_dB    = frequencyData{:,2};   % dB
phase_deg = frequencyData{:,3};   % degrees

% Convert to complex response
mag_lin  = 10.^(mag_dB/20);            % dB -> linear
phase_rad = deg2rad(phase_deg);       % degrees -> radians
resp = mag_lin .* exp(1i*phase_rad);  % complex response (Nx1 for SISO)

% FRD expects frequency in rad/s by default, so convert Hz -> rad/s
w_rad = 2*pi * f_Hz;

% Create frd (continuous-time SISO)
sys = frd(resp, w_rad);

% Record that the original frequency units were Hz
sys.FrequencyUnit = 'Hz';


% Details about Estimation Data
 % Import   sys

% State space model estimation
 Options = ssestOptions;
 Options.Display = 'off';  % Changed from 'on', which caused an error here on Answers
 Options.Focus = 'simulation';
 Options.InitialState = 'zero';
 Options.SearchOptions.Tolerance = 0;
 Options.SearchOptions.MaxIterations = 0;
 Options.N4Horizon = [15 15 15];

 ss1 = ssest(sys, 7, 'Form', 'free', Options);

 % bode plot
 % Get frequency vector from frd (rad/s)
w = 2*pi * sys.Frequency;   % if sys.Frequency is in Hz

opts = bodeoptions;
opts.PhaseWrapping = 'on';

% Plot both on the same Bode plot, specifying line styles
figure
% bode(ss1, w, sys, opts);
bodeplot(ss1,sys,opts)
grid on
legend('ss1','measured frd');

% Finally do step response
step(ss1);

[Arkadiy Turevskiy]

Hi Paul, I am not sure why the bode plot is not working for you. I just reuploaded my code to make sure I did not mess up somewhere, but it runs fine for me.

Both bode and bodeplot work just fine:

 % bode plot
 % Get frequency vector from frd (rad/s)
w = 2*pi * sys.Frequency;   % if sys.Frequency is in Hz

opts = bodeoptions;
opts.PhaseWrapping = 'on';

% Plot both on the same Bode plot, specifying line styles
bode(ss1, w, sys, opts);
grid on
legend('ss1','measured frd');

% Plot both on the same Bode plot, specifying line styles
figure
% bode(ss1, w, sys, opts);
bodeplot(ss1,sys,opts)
grid on
legend('ss1','measured frd');

I do see that when I try to run the code in the answer itself it errors out - I will need to look into that. It runs fine on the desktop MATLAB.

[Paul]

I was able to get the code to run on Answers as shown above by changing the display option in this line

Options.Display = 'off';  % Changed from 'on', which caused an error here on Answers

The Bode plot here on Answers is .... bizarre. Thanks for posting the pictures from your desktop run.

I don't undserstand this line:

bode(ss1, w, sys, opts);

Why would w be an input argument? In that syntax, wouldn't w be treated as an LTI object? But w is not displayed on the plot, so I guess I'm confused about that syntax.

I think this is an interesting approach and a nice match is obtained over the frequency range for which the data is provided. But the identified model also shows a deep null at the higher frequencies where there is no data, and we have no idea if that feature is justified or not.

[Arkadiy Turevskiy]

Thanks Paul for your comments and for figuring out the issue was in that line. I will check with Answers team about this behavior.

For the bode call that is confusing you - I had MATLAB Copilot create it for me. After your question I went through the doc and I agree it is not very clear that the order of arguments like in the call would work. So I asked MATLAB Copilot to explain and here is what it said.

Finally, to your third point - I agree that we have no idea whether high frequency behavior of the model is justified or not. I think ideally you woud have some idea about the physics of the system - what order dynamics it has. You can then select that model order when you estimate the data to have more physically meaningful results.

In this case as we don't know anything about what sort of system this data is for we can try estimating models of different orders.

I tried 2nd order and it was not a good fit at all.

3rd order is starting to look better.

See below for bode plot comparing 3rd order system - sys3, 7th order system - sys7. I also included step response plot.

[Paul]

Now that you (or actually Matlab Copilot) mention it, I think I've run into this before where a Control Systems Toolbox function was clever enough to disambiguate the arguments that were out of order. Maybe it was lsim.

But it is a bit surprising (to put it nicely) that Matlab Copilot would return code using an undocumented feature, when it could just as well return the code with the arguments in the documented order, as it actually shows as a "Recommendation," which MC doesn't follow.

Wholeheartedly agree that it's best to have an idea of the underlying physics/dynamics of the system to inform the system identification process.

[Arkadiy Turevskiy]

Another good point. MATLAB Copilot, like any AI tool, should be used carefully. I went back to the chat panel and flagged this recommendation so that it can rever to recommending the documented way to use the function. Thanks

[Andrew Ouellette]

@Paul If you're curious, the parser is written in plain m code:

str = which("parseRespFcnInputs");
type(str)

function [sysList,ParamList,PlotOptions] = ...
    parseRespFcnInputs(InputList,InputNames)
% Parser for input argument list to response plot function.

%   Copyright 1986-2024 The MathWorks, Inc.

PlotOptions = [];
ParamList = cell(1,0);
% System = system object, Name = system name, Style = plot style
sysList = struct('System',cell(0,1),'Name',[],'Style',[]);

% Scan input list
ni = nargin;
nsys = 0;
for ct=1:length(InputList)
    arg = InputList{ct};
    if isa(arg,'DynamicSystem')
        if isa(arg,'idlti') && size(arg,2) == 0
            arg = noise2meas(arg);
        end
        nsys = nsys+1;
        sysList(nsys).System = arg;
        if isempty(arg.Name) && ni>1
            sysList(nsys).Name = InputNames{ct};
        else
            sysList(nsys).Name = arg.Name;
        end
    elseif isa(arg,'plotopts.PlotOptions')
        PlotOptions = arg;
    elseif (ischar(arg) || isStringScalar(arg)) && ~any(strcmpi(arg,{'inv','zoh','foh','x0'}))
        % Note: support step(sys1,sys2,'r',sys3,sys4,'g--')
        if nsys==0 || ~isempty(sysList(nsys).Style)
            % Plot style comes first or is repeated
            error(message('Control:analysis:rfinputs01'))
        end
        % Validate plot style
        [~,~,~,msg] = colstyle(arg);
        if ~isempty(msg)
            error(message('Control:analysis:rfinputs02',arg))
        end
        sysList(nsys).Style = arg;
    elseif isa(arg,'iddata')
        % IDDATA no longer supported
  	  error(message('Controllib:plots:IddataNotSupported'))
    else
        ParamList = [ParamList {arg}]; %#ok<*AGROW>
    end
end

Basically, we extract all (model,linespec) pairs from the input list and pass all other arguments into a second parser later. As long as the remaining inputs are in the right order, you won't receive any errors. That being said, I wouldn't rely on this behavior since the parser could become stricter in the future.