calculate frequency response from two signal on domain time in matlab

Hi, I have arrays sample in time domain of input and output signal. size of vector is 2002 samples. I want to find frequency response. I wrote [h,w] = freqz(y,x,2002); figure() plot(abs(h)) but I get very strange graph. I assume to see filter. I want to calculate frequency response of smartphone(Audio filter) What do I do wrong ?

[Birdman]

Can you share your file?

[Alex Dashevsky]

x is input and y is output

[h,w] = freqz(x,y,2002); plot(abs(h))

Answer 1

[Star Strider]

You cannot get any useful information from your data because of the noise. However, with the correct sampling frequency (you did not state that or give a time vector), this will give you all the information it is possible to get from your data:

D1 = load(x.mat);
D2 = load(y.mat);
x = D1.x;
y = D2.y;

L = numel(x);                                   % Length Of Signal
Fs = 1;                                         % Sampling Frequency
Fn = Fs/2;                                      % Nyquist Frequency
FTx = fft(x)/L;
FTy = fft(y)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;             % Frequency Vector
Iv = 1:length(Fv);                              % Index Vector

TFyx = FTy./FTx;

figure(1)
subplot(2,1,1)
semilogx(Fv, 20*log10(abs(TFyx(Iv))))
ylabel('|H(f)| (dB)')
subplot(2,1,2)
semilogx(Fv, angle(TFyx(Iv))*180/pi)
ylabel('Phase (°)')
xlabel('Frequency (Hz)')

Note that you are only using one frequency, so if you want the frequency response at that frequency, divide the amplitude of ‘y’ by the amplitude of ‘x’. If you want the frequency response, the input must either be a sweep signal (such as a chirp) or broadband noise (essentially an impulse response). When you do that, my code will work for a broadband input signal producing a broadband output signal.

[Alex Dashevsky]

Could you explain with example ? The graph seem as noise, not filter. I enter sinus to my phone. I want to calculate H(f).

[Star Strider]

With only one frequency, you can calculate the transfer function only at that frequency. Simply divide the amplitude of the output by the amplitude of the input to get |H(f)| at that frequency.

If you want the broadband frequency response, for example from 20 Hz to 20000 Hz, you have to provide an input signal that goes from 20 Hz to 20000 Hz, and record the output at each frequency.

A broadband noise signal is the easiest way to do this. You have to use the output for the same signal you use as input, or the result will not be valid.

[Alex Dashevsky]

You divide it in your code. Right ? I have echo application. The output is after echo. The application sample with 48Khz. The problem is that I assume to see filter , not noise.

[Star Strider]

Yes. I calculate the complex transfer function as:

TFyx = FTy./FTx;

When I plotted your data in the time domain, there were two signals of different amplitudes that were 180° out of phase.

If you want to do echo analysis, the cceps (link), rceps (link) or related functions might be what you want. However, when I did that analysis, I saw no echo:

L = numel(x);                                   % Length Of Signal
Fs = 4.8E+4;                                    % Sampling Frequency
Fn = Fs/2;                                      % Nyquist Frequency
t = linspace(0, L, L)/Fs;                       % Time Vector
Ts = mean(diff(t));                             % Sampling Interval
FTx = fft(x)/L;
FTy = fft(y)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;             % Frequency Vector
Iv = 1:length(Fv);                              % Index Vector

TFyx = FTy./FTx;

figure(1)
subplot(2,1,1)
semilogx(Fv, 20*log10(abs(TFyx(Iv))))
ylabel('|H(f)| (dB)')
subplot(2,1,2)
semilogx(Fv, angle(TFyx(Iv))*180/pi)
ylabel('Phase (°)')
xlabel('Frequency (Hz)')

figure(2)
subplot(2,1,1)
semilogx(Fv, 20*log10(abs(FTx(Iv))))
ylabel('\itF\rm\{x(t)\}')
grid
subplot(2,1,2)
semilogx(Fv, 20*log10(abs(FTy(Iv))))
ylabel('\itF\rm\{y(t)\}')
xlabel('Frequency (Hz)')
grid

figure(3)
plot(t, x,    t, y)
grid

rc = rceps(y);
cc = cceps(y);

figure(4)
subplot(2,1,1)
plot(t, rc)
axis([min(t)  max(t)    -0.1  +0.2])
title('Real Cepstrum')
subplot(2,1,2)
plot(t, cc)
axis([min(t)  max(t)    -2  +2])
title('Complex Cepstrum')

[Alex Dashevsky]

I only want to calculate H from input and output. I enter the sinus to my phone. the sinus enter to my audio application and play output signal.

[Star Strider]

With one input frequency you will (and in a linear system have gotten) one output frequency. The transfer function is therefore the Fourier transform of the output divided by the Fourier transform of the input. You have to do this at more than one frequency to get a transfer function over those frequencies.

[Alex Dashevsky]

Could you give example ?

[Star Strider]

I already did. See my Answer and Comments.

[Alex Dashevsky]

I don't understand. If I do this procedure on two frequency the I will get two H. What after ?

[Star Strider]

No. You will get one H(f).

This line calculates the complex transfer function:

TFyx = FTy./FTx;

and these lines plot the transfer function magnitude as ‘H(f)’ in the upper subplot, and the phase in the lower subplot:

figure(1)
subplot(2,1,1)
semilogx(Fv, 20*log10(abs(TFyx(Iv))))
ylabel('|H(f)| (dB)')
subplot(2,1,2)
semilogx(Fv, angle(TFyx(Iv))*180/pi)
ylabel('Phase (°)')
xlabel('Frequency (Hz)')

The complex transfer function is defined as the Fourier transform of the output (that I assume here is ‘y’) divided by the Fourier transform of the input (that I assume is ‘x’). You can then compute the magnitude (using the abs function) and phase (using the angle function) from the complex transfer function, as I have done here in my code.

If you want to calculate the frequency response of a system, you must present a wide range of frequencies to the input, and record the input to and output of the system at the same time. Then you can compute the transfer function from those time-domain data by doing the Fourier transform on both.

[Alex Dashevsky]

Sorry, but I can't understand how I will get only one H ? I on every frequency get vector of input samples and output samples. I use two kind of frequency. 1. frequency of input signal. 2 sample frequency. Do you mean 1 ? Can you give me example ?

[Star Strider]

I have already explained this.

‘H’ is the transfer function. There is only one, defined at every frequency, calculated by these equations:

FTx = fft(x)/L;
FTy = fft(y)/L;
TFyx = FTy./FTx;
H = abs(TFyx(Iv));
Phase = angle(TFyx(Iv))

If this does not answer your question, you need to be more specific in describing what you want.

[Alex Dashevsky]

You wrote that if I want to find frequency response of the system then I need to do this on range of frequency. Can you give example on two frequency how I find H ?

[Star Strider]

You are doing this in hardware (your cell phone), so I cannot demonstrate it.

If you want to generate an input signal, use this:

t = linspace(0, 4.8, 4.8E+4);
f0 =    20;
f1 = 10000;
Fs = 1/mean(diff(t));

x = chirp(t,f0,t(end),f1);

sound(x, Fs)

Use ‘x’ here as your input and your cell phone recording of it as your output, ‘y’.

That is the best I can do.

[Omar]

Hello Star Strider,

Could you please advise me why you chose Nyquist Frequency (Fn) to clacualte the frequency vector not the sampling frequency Fs itself in this example.

Please tell me.

[Nghi Nguyen]

Hello Star Strider,

May I have a question about the unit of the x-axis? Why is it in Hz, but not in rad/s?

Thank you very much.

Nghi