CODEWITHSUNDEEP
matlabsignal-processingfftconvolution
yankeefan11
yankeefan11

Reputation: 485

Fourier Transforms in MatLab

So I have had a few posts the last few days about using MatLab to perform a convolution (see here). But I am having issues and just want to try and use the convolution property of Fourier Transforms. I have the code below: 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')

My goal is to take my case, which is $measured = rt \ast signal$ and find signal. Once I find my signal, I convolve it with the rectangle and should get back measured, but I do not get that.

I have very little matlab experience, and pretty much 0 signal processing experience (working with DFTs). So any advice on how to do this would be greatly appreciated!

Upvotes: 0

Views: 1107

Answers (1)

macduff
macduff

Reputation: 4685

After considering the problem statement and woodchips' advice, I think we can get closer to a solution.

Input: u(t)

Output: y(t)

If we assume the system is causal and linear we would need to shift the rect function to occur before the response, like so:

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

Shifted Input

Next, consider the response to the input. It looks to me to be first order. If we assume as such, we will have a system transfer function in the frequency domain of the form:

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

You had been trying to use convolution to and FFT's to find the impulse response, "transfer function" in the time domain. However, the FFT of the rect, being a sinc has a zero crossing periodically. These zero points make using the FFT to identify the system extremely difficult. Due to:

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

So we have U(s) = A*sinc(a*s), with zeros, which makes the division go to infinity, which doesn't make sense for a real system.

Instead, let's attempt to fit coefficients to the frequency domain linear transfer function that we postulate is of order 1 since there are no overshoots, etc, 1st order is a reasonable place to start.

EDIT I realized my first answer here had a unstable system description, sorry! The solution to the ODE is very stiff due to the rect function, so we need to crank down the maximum time step and use a stiff solver. However, this is still a tough problem to solve this way, a more analytical approach may be more tractable.

We use fminsearch to find the continuous time transfer function coefficients like:

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

Something like that anyway should get you going.

x = findTf([1 1 1]',rt',measured',x');

Upvotes: 2

Related Questions