Using IFFT for obtaining time response of measured freq response

18 次查看(过去 30 天)
Hello, I have numerical frequency response data (G(s=j2*pi*f), f) for a system. I want to obtain the impulse response in time domain. I should be able to obtain this using the IFFT function. But the scaling is not clear to me. Here is the code I am using. What am I missing? Thanks. Amit.
clc; close all;
fmax = 1e6;
L = 1024;
fdelta = 2*fmax/(L-1)
fgrid = [-fmax:fdelta:fmax]; %fmax = fs/2 => deltaT = 1/(2*fmax)
tau = 1/(2*pi*5e4)
% Y(s) = 1/(1+s.tau) Given Freq Domain data
Y = 1./(1+j*2*pi*fgrid*tau);
y_t = ifft((Y));
Tdelta = 1/(2*fmax);
t = [0:Tdelta:Tdelta*(length(y_t)-1)];
% compare with known time domain function
figure; plot(t, y_t, 'bx', t, exp(-t/tau), 'r');
grid; axis([0 1e-4 0 1])

回答(4 个)

Rick Rosson
Rick Rosson 2011-12-28
Please try:
Fs = 2e6;
L = 1024;
dF = Fs/L;
f = (-Fs/2:dF:Fs/2-dF)';
s = j*2*pi*f;
Fc = 50e3;
alpha = 2*pi*Fc;
tau = 1/alpha;
Y = alpha./(s+alpha);
dt = 1/Fs;
t = dt*(0:L-1)';
y = L*ifft(ifftshift(Y));
x = exp(-t/tau);
HTH.
Rick
  2 个评论
Amit
Amit 2011-12-28
Rick,
Thanks for the reply.
I tried the code you suggested. I still get a scaling error.
>> [(x(1:10)) abs(y(1:10))']
ans =
1.0000e+000 77.8669e+000
854.6360e-003 152.1368e+000
730.4027e-003 109.5818e+000
624.2284e-003 105.7690e+000
533.4881e-003 81.7647e+000
455.9381e-003 76.5834e+000
389.6611e-003 59.9684e+000
333.0184e-003 55.8890e+000
284.6095e-003 43.7460e+000
243.2376e-003 40.9324e+000
>> [(x(1:10))./abs(y(1:10))']
ans =
12.8424e-003
5.6176e-003
6.6654e-003
5.9018e-003
6.5247e-003
5.9535e-003
6.4978e-003
5.9586e-003
6.5060e-003
5.9424e-003
>> 2*pi
ans =
6.2832e+000
>> 2*pi/1024
ans =
6.1359e-003
Looks like there should be scaling factor of 2*pi instead of L=1024.
However, increasing the sampling frequency changes the scaling.
With fs=200e6.
>> [(x(1:10))./abs(y(1:10))']
828.7821e-003
463.9740e-003
517.4152e-003
484.3603e-003
507.4654e-003
489.3110e-003
504.0767e-003
491.5127e-003
502.3784e-003
492.7540e-003
I am really not sure what is happening. I wish there were a matlab function to take care of all this and just give me a time domain impulse response!
Amit.
LI QZH
LI QZH 2016-12-22
编辑:LI QZH 2016-12-23
Y = alpha./(s+alpha); i think this is not the right Laplace transform
y=1./(s+alpha) is the right form
am i right ? thanks

请先登录,再进行评论。


Rick Rosson
Rick Rosson 2011-12-28
Do you have access to either the Control Systems Toolbox or the Signal Processing Toolbox? If so, which one (or both)?

Rick Rosson
Rick Rosson 2011-12-29
Please check this related answer:
HTH.

Rick Rosson
Rick Rosson 2011-12-29
I modified the code I posted earlier to correct the scaling factor:
Fs = 2e6;
L = 1024;
dF = Fs/L;
f = (-Fs/2:dF:Fs/2-dF)';
s = j*2*pi*f;
Fc = 50e3;
alpha = 2*pi*Fc;
tau = 1/alpha;
Y = alpha./(s+alpha);
dt = 1/Fs;
t = dt*(0:L-1)';
G = 2*pi;
y = G*abs(ifft(ifftshift(Y)));
x = exp(-t/tau);
figure;
plot(t,x,t,y);
HTH.
Rick
  1 个评论
Greg Heath
Greg Heath 2011-12-30
FFTing and IFFTing is tricky business because you have to be perfectly clear about original assumptions. For example, time domain signals contain N measurements at a sampling frequency
Fs. If the signal is real the corresponding nonnegative frequency spectrum measurements have a spacing df = Fs/N and
have either the range
f = df*[0:N/2], Neven
fmax = Fs/2
L = length(f) = N/2+1
with L odd when N is even
or
f = df*[0:(N-1)/2], N odd
fmax = Fs/2 - df/2
L = (N+1)/2
with L even when N is odd
Now since you seem to have L = 1024 measurements, my conclusion is that
N = 2*L + 1 = 2049 is odd
and
Fs = 2*N*fmax/(N-1) = 2.001e6
Hope this helps.
Greg

请先登录,再进行评论。

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by