0

所以最近几天我有一些关于使用 MatLab 执行卷积的帖子（见这里）。但我遇到了问题，只想尝试使用傅里叶变换的卷积属性。我有下面的代码：width = 83.66;

x = linspace(-400,400,1000);

      a2 =  1.205e+004  ;
       al =  1.778e+005  ;
       b1 =       94.88  ;
       c1 =       224.3  ;
       d =       4.077  ;

measured =  al*exp(-((abs((x-b1)./c1).^d)))+a2;

%slit
rect = @(x) 0.5*(sign(x+0.5) - sign(x-0.5));
rt = rect(x/width);


subplot(5,1,1);plot(x,measured);title('imported data-super gaussian')
subplot(5,1,2);plot(x,(real(fftshift(fft(rt)))));title('transformed slit')
subplot(5,1,3);plot(x,rt);title('slit')

u = (fftshift(fft(measured)));
l = u./(real(fftshift(fft(rt))));
response = (fftshift(ifft(l)));

subplot(5,1,4);plot(x,real(response));title('response')

%Data Check
check = conv(rt,response,'full');
z = linspace(min(x),max(x),length(check));
subplot(5,1,5);plot(z,real(check));title('check')

我的目标是接受我的案例，即 $measured = rt \ast signal$ 并找到信号。一旦我找到我的信号，我将它与矩形进行卷积并且应该返回measured，但我不明白。

我的 matlab 经验很少，信号处理经验几乎为零（使用 DFT）。因此，任何有关如何执行此操作的建议将不胜感激！

4

1 回答 1

2

在考虑了问题陈述和 woodchips 的建议后，我认为我们可以更接近解决方案。

输入：u(t)

输出：y(t)

如果我们假设系统是因果和线性的，我们需要将rect函数转换为在响应之前发生，如下所示：

rt = rect(((x+270+(83.66/2))/83.66));
figure; plot( x, measured, x, max(measured)*rt )

移位输入

接下来，考虑对输入的响应。在我看来，这是第一顺序。如果我们这样假设，我们将在以下形式的频域中具有系统传递函数：

H(s) = (b1*s + b0)/(s + a0)

您一直在尝试使用卷积和 FFT 来找到脉冲响应，即时域中的“传递函数”。但是， 的 FFTrect是 asinc周期性地过零。这些零点使得使用 FFT 识别系统变得极其困难。由于：

Y(s)/U(s) = H(s)

所以我们有U(s) = A*sinc(a*s), 带零，这使得除法变为无穷大，这对于实际系统没有意义。

相反，让我们尝试将系数拟合到我们假设为 1 阶的频域线性传递函数，因为没有过冲等，1 阶是一个合理的起点。

编辑 我意识到我在这里的第一个答案有一个不稳定的系统描述，对不起！由于rect函数的原因，ODE 的求解非常僵硬，因此我们需要降低最大时间步长并使用刚性求解器。然而，这仍然是一个很难解决的问题，一个更具分析性的方法可能更容易处理。

我们fminsearch用来找到连续时间传递函数系数，例如：

function x = findTf(c0,u,y,t)
  % minimize the error for the estimated
  % parameters of the transfer function
  % use a scaled version without an offset for the response, the
  % scalars can be added back later without breaking the solution.
  yo = (y - min(y))/max(y);
  x = fminsearch(@(c) simSystem(c,u,y,t),c0);
end

% calculate the derivatives of the transfer function
% inputs and outputs using the estimated coefficient
% vector c
function out = simSystem(c,u,y,t)
  % estimate the derivative of the input
  du = diff([0; u])./diff([0; t]);
  % estimate the second derivative of the input
  d2u = diff([0; du])./diff([0; t]);
  % find the output of the system, corresponds to measured
  opt = odeset('MaxStep',mean(diff(t))/100);
  [~,yp] = ode15s(@(tt,yy) odeFun(tt,yy,c,du,d2u,t),t,[y(1) u(1) 0],opt);
  % find the error between the actual measured output and the output
  % from the system with the estimated coefficients
  out = sum((yp(:,1) - y).^2);
end

function dy = odeFun(t,y,c,du,d2u,tx)
  dy = [c(1)*y(3)+c(2)*y(2)-c(3)*y(1); 
        interp1(tx,du,t);
        interp1(tx,d2u,t)];
end

无论如何，这样的事情应该让你继续前进。

x = findTf([1 1 1]',rt',measured',x');
于 2013-07-02T15:11:30.783 回答