MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn moreOpportunities for recent engineering grads.

Apply Today
Asked by buscuit on 15 May 2012

Hi to everybody,

I have implemented a little function that calculates the impulse response from a complex dataset (FR_data).

I am trying to test how different FR_data influnce the impulse response so I plot the FR_data over frequency using plot(frequency,abs(FR_data)) and after I ran the routine and plot the results again using plot(time,abs(y)).

the routine looks like:

function y=IR(FR_data)

L=length(FR_data); Y=zeros(1,2*length(FR_data)); Y(2:L+1)=FR_data; Y(L+2:end)=fliplr(FR_data(1:end-1)); y=ifft(Y); end

I notice that when I boost the higher frequencies in my FR_data the peak of the impulse response gets closer to time 0 and increases in amplitude and is somehow shorter. When I boost the low frequencies then the peak of the impulse response gets further away from time 0, decreases in amplitude and is wider.

Can anybody enlighten me please on if this is normal behaviour?

Furthermore, all the impulse responses exhibit a very very short and sharp maximum amplitude on time 0 (actually it extends for 2-3 samples, then dies out and then the impulse response follows. Is there a way to interpret this?

Any input would be greatly appreciated.

thank you for your patience

A small addition that might explain better what I am talking about. My datasets looks like: https://picasaweb.google.com/106674205971471180972/Matlab?authkey=Gv1sRgCKjU47ij4922SQ#5743069984065489538 they represent complex data that extend from 1 Hz to 20 kHz with a step of 1 Hz. So 20000 data. When I plot abs of those in a semilogx I get: https://picasaweb.google.com/106674205971471180972/Matlab?authkey=Gv1sRgCKjU47ij4922SQ#5743069982822662066. They represent 3 different FR_data sets. The x axis represents frequency.

After I run my routine I get IRs that look like: https://picasaweb.google.com/106674205971471180972/Matlab?authkey=Gv1sRgCKjU47ij4922SQ#5743069982995627554 and a zoomed version in the beginning of the x axis to actually see the IRs: https://picasaweb.google.com/106674205971471180972/Matlab?authkey=Gv1sRgCKjU47ij4922SQ#5743070000860043586 The x axis here represents the samples (not the time).

*No products are associated with this question.*

Answer by buscuit on 15 May 2012

Hi Elige,

thanks for the contribution once again. I would like to upload some images on what i am talking about here, indeed, but I can't seem to find a way to attach them on my post.

Elige Grant on 15 May 2012

You can do this by uploading a picture to the web (I upload to my Picasa website) and then you just put the picture URL inside << www.something.com/example/mypic.jpg >>. You can add this as an edit to the bottom your original question. Let me know if that doesn't work.

Elige Grant on 16 May 2012

Those impulse responses do look a bit strange. Are those still plotted using "plot(time,abs(y))" or "plot(time,y)" with no ABS? I think it will be easier for me to diagnose if you either A) copy/paste the complex FR_data into another answer here or B) supply the function that was used to create the FR_data.

Answer by buscuit on 16 May 2012

The IR plots are plotted using plot(abs(y)). So no time on x axis just samples.

The function I am using to generate my FR_data is:

function [G,f]=planeFR(wx,wy,h,res)

c=343; f=1:res:20000; B=1; R12=1; R1=h; R2=h; xl=wx/2; yl=wy/2; k=2*pi*f/c; psi=0; l=1./k;

% syms u v;

% Pinf=zeros(1,length(k)); % P=zeros(1,length(k));

% for i=1:1:length(k)

ul= cos(psi) * sqrt( (R1+R2)./(2*l*R1*R2) ).* wx; vl= sqrt( (R1+R2)./(2*l*R1*R2) ).* wy;

% A=(1j*B*R12*exp(-1j*k(i)*(R1+R2)))/(2*(R1+R2)); % I1=int(exp(-1j*(pi/2)*(u^2)),u,-ul,ul); % I2=int(exp(-1j*(pi/2)*(v^2)),v,-vl,vl);

[Cul Sul]=fcs(ul); [Cvl Svl]=fcs(vl);

I1=2*( Cul - 1j*Sul); I2=2*( Cvl - 1j*Svl);

G = ( 1j* I1.* I2 )/2;

% semilogx(k,abs(G),'r');

end

where fcs is taken from http://www.mathworks.com/matlabcentral/fileexchange/6580 in order to calculate the fresnel integrals.

I use h=2, res=1 and I play with wx and wy.

What this function gives me is the frequency response I was talking about.

Elige Grant on 16 May 2012

Can you give me an example of the "wx" and "wy" you used to create those figures... just so I am in the same ballpark range.

Answer by Elige Grant on 17 May 2012

I reran through the process of generating the Frequency Response data (G) using your code and the one from the FEX just to see if there could have been any steps that were not properly handled in converting it to the time-domain... but I still got the same bizarre looking Impulse Response you provided pictures of. As I am not familiar with the type of FR function you are using, I cannot really comment on why the shape looks so bizarre... but I can tell you that the "spike" at your first time-domain data sample is basically equal to **mean(real(G))**. This is simply a consequence of having real frequency amplitudes that have a static shift, in this case they are all pretty much centered around 1... so you get a spike with amplitude 1 at your first time-domain data sample. If you were to subtract your frequency-domain data by 1, you would get rid of this spike in the time-domain.

Not sure if I have any other helpful advice if you are expecting a different shaped IR in the time domain without knowing where this FR formulation is coming from.

***********************

**EDIT**

**********************

I forgot to mention, that I would change your :

f=1:res:20000;

to:

f=res:res:20000;

because your first non-zero frequency will not always be 1 (only when res=1 is that valid).

I am only going to add the next bit just so you can double check the way you create your symmetric frequency response data (a little more general than my previous post):

If:

[G,f]=planeFR(wx,wy,h,res);

Then:

N = numel(G); % number of non-zero, positive frequency amplitudes df = (f(end)-f(1))/(N-1); % frequency increment Nyq = f(end); % Nyquist frequency Fs = 2*Nyq; % Sample frequency %Fs = 2*N*df; %<--- should be same as above newN = 2*N; % Number of samples in the full frequency response spectrum %newN = Fs/df; %<--- should be same as above

newG = zeros(newN,1); % Pre-allocate full spectrum array newG(2:1+N) = G; % populate non-zero, positive frequency amplitudes newG(2+N:end) = fliplr(conj(G(1:end-1))); % neg. frequency amplitudes % newG = newG - mean(real(newG)); %<---- remove that spike in time-domain % newG = newG*Fs; %<---- correct time-domain amplitudes newf = -Nyq : df : Nyq-df;

% Plot Frequency Response % figure; plot(newf, abs(newG));

% Plot Impulse response t = (0:newN-1)/Fs; % time g = ifft(newG); figure; plot(t,g);

buscuit on 17 May 2012

Elige, first of all, I would like to thank you for all the time you spent on my problem.

My FR_data give me the frequency response of the reflections from a surface when hit by a sound wave. They are normalized around 1, that is why all the frequency responses converge to that point.

When the surface is big (wx,wy) compared to the distance (h) of the source of the sound wave from the object the frequency response gets its maximum at low frequencies. When wx, wy are small compared to h then the frequency response is shifted towards higher frequencies.

I was asked to get the impulse responses out of those frequency responses, but I was not able to interpret them. I think I understand what you said about the spike at the beginning of the IRs but still, I cannot understand why the second spikes are getting further and further from the origin as wx,wy increasing (or in other words when the FRs are extending on a wider range of frequency).

Do you have any recommendation on any reference I could study in order to understand more on how to interpret IRs?

## 1 Comment

Direct link to this comment:http://www.mathworks.es/matlabcentral/answers/38465#comment_79565

A note of clarification (only because I helped in a previous post: http://www.mathworks.com/matlabcentral/answers/38316-obtain-the-impulse-response-from-a-frequency-response) is that FR_data represents the non-zero, positive frequency amplitude values in the frequency domain.

I would make one small modification to your function (just so it is more applicable to other situations in the future). Change your line from "fliplr(FR_data(1:end-1))" to "fliplr(conj(FR_data(1:end-1)))" in order to account for the possibility that FR_data is complex (instead of only real valued).

The fact that your pulse gets narrower with a larger contribution of high frequencies or wider with a larger contribution of low frequencies seems normal. But.. is it possible that you post a picture of "plot(time,y)" (without the "abs") or copy/paste the values associated with a high-frequency and low-frequency version of "FR_data" ?