Need some matlab help!! / Computer Science / Forums

Forums

4hv.org :: Forums :: Computer Science

« Previous topic | Next topic »

Need some matlab help!!

1 2 next

Move Thread

LAN_403

scrooch

Mon May 10 2010, 08:37AM

Registered Member #908 Joined: Wed Jul 18 2007, 05:53AM
Location:
Posts: 49

I have to implement
G(v)-W(v)

where G(v) and W(v) are in the time domain(omega instead of v would have been more appropriate)

I have the samples for G(v) in time domain in a variable called Gt
so to convert it to frequency domain i used the fft function in matlab as
Gv=fft(Gt)

i do the same for W(v)

now if i go
Gv-Wv in matlab is it the same as G(v)-W(v) in the mathematical sense?
as in do i have to mess around with the fft function or is it as simple as above?

thanks for your input:-)

Dr. Slack

Mon May 10 2010, 06:05PM

Registered Member #72 Joined: Thu Feb 09 2006, 08:29AM
Location: UK St. Albans
Posts: 1659

I'm not clear that I understand the assertion that G(v) and W(v) are in the time domain.

Is it the case that you have a function G, which you have got specifed in the time domain as a vector of samples Gt. You can also specify it in the frequency domain as a function of v as G(v), which you could call Gv in matlab. Getting from the time domain to the frequency domain may be as simple as performing an FFT of the time samples vector to get a frequency samples vector, or it may be more complicated, depending on what you mean by G(v). Any difficulties you may be having may be to do with the time to frequency process, or they may be something else.

Matlab's Gv-Wv will do exactly the same as G(v) - W(v) in the mathematical sense. However, is the subtraction of two vectors of complex numbers that you want? To help answer this question, what is f(v)? Is it ...

a) the complex voltage spectrum of an unwindowed time vector (what you get out of a raw FFT)
b) the power spectrum of a windowed time vector (the sort of spectrum a spectrum analyser produces)
c) or any one of a number of alternative spectra definitions

The answer may be obvious within any particular field of application - deconvolution, frequency response normalisation, or answering homework assignment question would all imply different approaches. What are you trying to do?

I don't know whether this helps, it's a little bit of matlab code demonstrating that the FFT is a linear function, in the sense that f(a-b) = f(a)-f(b)

%% demonstrates fft linearity

clear;
clc;

% make random time vectors
Gt=randn(10,1);
Wt=randn(10,1);
% take their FFTs
Gv=fft(Gt);
Wv=fft(Wt);
%difference the FFts
Dv1=Gv-Wv;
% alternatively fft their difference
Dv2=fft(Gt-Wt);
% is there any difference in the two routes?
display(Dv1-Dv2);

and the last line prints out the answer

ans =

  1.0e-015 *

        0          
        0 + 0.2220i
        0 + 0.4441i
  -0.2220 - 0.4441i
  -0.4441 - 0.2220i
   0.4441          
  -0.4441 + 0.2220i
  -0.2220 + 0.4441i
        0 - 0.4441i
        0 - 0.2220i

so the difference is at the level of double precision rounding noise. QED, but I'm not sure you wanted that demonstrated.

scrooch

Tue May 11 2010, 07:55AM

Registered Member #908 Joined: Wed Jul 18 2007, 05:53AM
Location:
Posts: 49

Hey Dr slack thanks for your quick reply. I actually had a typo in my post
W(v) and G(v) are in the frequency domain and W(t) and G(t) are in time domain sorry about that. So i transformed them from time to frequency domain using fft function.

This is for a blind deconvolution problem

I will restate my question again in-terms of the actual variables.

g(t)= f(t)*d(t) + w(t) //in time domain

G(v)=F(v)D(v) + W(v) //in frequency domain

where g(t) is the original signal that has to be cleared up using deconvolution.
f(t) is the actual data.
w(t) is noise
d(t) is the impulse function
* denotes convolution
the gist of my question is( i want an idea of how the fft function works)

can
D(v)= |G(v) - W(v)| //(the actual equation to be implemented has some constants as well which i did not include here for simplicity.)

be implemented in matlab as

D(v) =|fft(Gt) - fft(Wt)| ?

Assuming i can do so i attempted it and got some samples for D(v) 64000 to be exact(The original signals were sampled at 16000Hz).

How can i visualize this transfer function in matlab(to see its stability, amplitude, phase response etc)?

Dr. Slack

Tue May 11 2010, 10:15AM

Registered Member #72 Joined: Thu Feb 09 2006, 08:29AM
Location: UK St. Albans
Posts: 1659

How can i visualize this transfer function in matlab(to see its stability, amplitude, phase response etc)?

you can visualise it using plot(func). Matlab can do all sorts of plots, I've only explored a tiny fraction of them. Type "help plot"

the gist of my question is( i want an idea of how the fft function works)

In that case, it's probably not best to start with the process of blind deconvolution. The easiest way to get a feel for what the FFT does is is to generate some test time vectors, FFT them, and plot the results. Make sure you understand what the plot routine is doing with the complex numbers the fft produces.

g(t)= f(t)*d(t) + w(t) //in time domain

G(v)=F(v)(?)D(v) + W(v) //in frequency domain

This is the sort of notation in text books that doesn't just drop into Matlab speak, I presume you're missing a way of typing the (?) convolution operator.

Unfortunately, there's quite a gulf between hearing about the fft function and using it for blind deconvoluton, you're trying to build a cabinet (blind deconvolution) before you've figured out how to use one of the tools, a hammer (fft).

a) You can use the fft function to implement convolution ( and deconvolution), but not like that. If you are in a matlab environment, it's fewer keystrokes, and cooler, to use the conv function.

if you type "help conv" you get, amongst other information, this

Algorithm
The convolution theorem says, roughly, that convolving two sequences
is the same as multiplying their Fourier transforms. In order to make this
precise, it is necessary to pad the two vectors with zeros and ignore roundoff
error. Thus, if X = fft([x zeros(1,length(y)-1)])and Y = fft([y zeros(1,length(x)-1)])then conv(x,y) = ifft(X.*Y)

the key word in that is "roughly"

there's also a deconv function

b) the fft function doesn't transform "time signals" to "frequency signals" and vice versa. It does transform from the time domain to the frequency domain and vice versa sure enough, but it transfroms a sampled cyclic aliased signal from one to the other. You have to take care of all of the issues of sampling, de-aliasing and making the signals cyclic yourself, outside the transform itself. If you don't, it will assume they are anyway, with usually unexpected and signal destroying results.

There's a lot to learn, the best way is by experiment, and Matlab's a fairly good vehicle for creation and visualisation of tests. Do find out what's involved in the theory of blind deconvolution though, another tool a cabinet maker needs is a plan of the cabinet. Good luck.

scrooch

Wed May 12 2010, 11:33AM

Registered Member #908 Joined: Wed Jul 18 2007, 05:53AM
Location:
Posts: 49

Thanks for the quick response again Dr Slack.
I read through the conv function and I did the padding of the ffts before i multiplied them.
I tested this(like in the matlab help page) by taking two sequences and then convolving them, and then again using fft.
They gave the same result and now ive got that bit clear.
I have one final issue. This was a blind deconvolution problem.

Now I have estimated the transfer function of the unknown filter(I have it's discrete values in frequency space).
What is the best method of deconvolving the degraded signal with this information?
I tried looking into the wiener filter but other than lots of theory i couldn't find any information on how to put it into practice. Are there other better suited methods?

Steve Conner

Wed May 12 2010, 12:18PM

Registered Member #30 Joined: Fri Feb 03 2006, 10:52AM
Location: Glasgow, Scotland
Posts: 6706

Blind deconvolution is a difficult problem. There's no general solution that I know of.

Techniques for doing it rely on assumptions about the nature of the signal you're trying to retrieve, or the transfer function you're trying to deconvolve, or both. Essentially, you make a model of the unwanted transfer function using prior knowledge, and estimate the parameters of your model from the signal, by looking at the difference between what the signal actually is like, and what you expect it ought to be like.

An example of what I mean is adaptive echo cancellation, which assumes that the unwanted transfer function is just a delay. And the LMS algorithm, which assumes that the signal you're interested in is random noise.

scrooch

Wed May 12 2010, 12:31PM

Registered Member #908 Joined: Wed Jul 18 2007, 05:53AM
Location:
Posts: 49

Steve McConner wrote ...

Blind deconvolution is a difficult problem. There's no general solution that I know of.

Techniques for doing it rely on assumptions about the nature of the signal you're trying to retrieve, or the transfer function you're trying to deconvolve, or both. Essentially, you make a model of the unwanted transfer function using prior knowledge, and estimate the parameters of your model from the signal, by looking at the difference between what the signal actually is like, and what you expect it ought to be like.

An example of what I mean is adaptive echo cancellation, which assumes that the unwanted transfer function is just a delay. And the LMS algorithm, which assumes that the signal you're interested in is random noise.

Thanks for the reply Steve. Im fiddling around with a blind deconvolution scheme.

From that scheme i have approximated the transfer function of the unknown filter the truth signal was sent through.

I tested my transfer function in a crude way.
Which was to do

G(v)/D(v)
where G(v) is the fourier transformed version of the dirty signal
And D(v) is the samples of the transfer function in frequency space.

The results were way better than what i thought. The resultant signal was way clearer.

But as you can see this "Inverse filtering" method is very crude.

So its this point where i need help. Is there any other better way where i can do the actual deconvolution?.

Steve Conner

Wed May 12 2010, 01:02PM

Registered Member #30 Joined: Fri Feb 03 2006, 10:52AM
Location: Glasgow, Scotland
Posts: 6706

The hard part of blind deconvolution isn't the deconvolution, that's mathematically very easy, even if it is computationally intensive. The method you mentioned (FFT the signal, divide by the transfer function, inverse FFT the result) works fine, provided you take care with windowing and so on.

The hard part is the "blind" - figuring out what the transfer function for deconvolution is.

Dr. Slack

Wed May 12 2010, 01:07PM

Registered Member #72 Joined: Thu Feb 09 2006, 08:29AM
Location: UK St. Albans
Posts: 1659

Unfortunately you are into PhD territory here, clever people spend years of their lives trying to improve the state of the art.

The two main problems (and I am not an expert of the subject) as I see it are ...

a) any noise contamination cannot be corrected, the whole point of noise is that it is unpredictable. A Weiner filter makes the best fist at estimating the most likely result. Any knowledge of characteristics of the signal can improve things.

b) if the corrupting system response has a zero or near zero in the frqeuency domain, then you cannot do the simple "divide by inverse of system frequency response" thang, as the corrected signal will go to infinity. Deciding what to do with the data lost down the zero depends on the situation. In medical imaging scenarios, you would just keep it zero so as not to add artifacts. In delivering entertainment content, you'd try to mask it in as imperceptible way as possible.

It all comes down to knowing how much you know about the underlying signal, how much you know about the interference, and what defects you can tolerate in the corrected results. It's perhaps more important to know when your algorithm has done as well as is theoretically possible, given the uncertainties, because then you can stop striving to improve it.

Steve Conner

Wed May 12 2010, 01:28PM

Registered Member #30 Joined: Fri Feb 03 2006, 10:52AM
Location: Glasgow, Scotland
Posts: 6706

Dr. Slack wrote ...

Deciding what to do with the data lost down the zero depends on the situation.

Or as it was put by the lecturer on a DSP course my work sent me to:
"Zero divided by zero equals anything you want."

If I remember right, "The Scientist and Engineer's Guide To DSP" contains a beginner's guide to this kind of stuff, and it's available online for free.

1 2 next

Moderator(s): Chris Russell, Noelle, Alex, Tesladownunder, Dave Marshall, Dave Billington, Bjørn, Steve Conner, Wolfram, Kizmo, Mads Barnkob