Kindly verify my code of DTMF decoding step

I want to decode DTMF tones by using FIR Filter. The filters that are required to be used in filter bank are cnstructed by Sinusoidal impulse response of the form h[n]= (2/L)cos(2*pi*fb*n/fs) where 0<n<L

L is filter length fb defines the freq location of passband e.g we can pick 697Hz

the book says to generate bandpass filter for 770 Hz component with L=50 and fs=12000. This has to be done by creating a vector of filter co-efficients ,h770 which are determined by evaluatiing above stated equation. plot the filter coefficients using stem().

I have done it in this way. Is it ok

h=[];
L=50;
fs=12000;
fb=770;
for n=1:L
    h(n)=(2/L)*cos(2*pi*fb*n/fs);
end
stem(h)

[Jan]

Please check the results with some test data.

Answer 1

[Wayne King]

Hi, Notice how Walter has indexed the for loop from 0, but was careful to index the vector h from 1.

h=[];
L=50;
fs=12000;
fb=770;
for n=0:L
    h(n+1)=(2/L)*cos(2*pi*fb*n/fs);
end
stem(h)

So that it is legal MATLAB indexing, however, I would say to avoid the for loop entirely. This is cleaner MATLAB.

L=50;
fs=12000;
fb=770;
% if you really want 0<n<L
h = (2/L)*cos(2*pi*fb*(1:L-1)/fs);
stem(h);

[Walter Roberson]

h = (2/L)*cos(2*pi*fb*(0:L-1)/fs);

for the 0 <= n < L case.

Answer 2

[moonman]

I dont have test data. i just want somebody to read my statement and see my code i hope it is ok but just want to cross check from experts

Answer 3

[Walter Roberson]

The description says 0 < n < L but you have n going from 1 to L which is one step too far (n == L)

[Walter Roberson]

Notice the strict inequality: 0 < n < L . That means that n cannot be 0 and cannot be L, so there is no need to use "for n=0:L" (not unless the instructions are wrong.)

If the instructions had said 0 <= n <= L, then you would code as

for n=0:L

h(n+1)=(2/L)*cos(2*pi*fb*n/fs);

end

[moonman]

The instructions say

0<=n < L (n is greater than or equal to 0)

How can i code this

???

[Walter Roberson]

for n=0:L-1

h(n+1)=(2/L)*cos(2*pi*fb*n/fs);

end

Answer 4

[moonman]

Yes, this issue i know but when i tell matlab to do

for n=0:L

It gives error

??? Attempted to access h(0); index must be a positive integer or logical.

How to remove this error

Answer 5

[moonman]

ok thanks king

kindly confirm me is my code ok as per book question

[Wayne King]

Yes, you can confirm by looking at the frequency response

fvtool(h,1,'Fs',fs);

Answer 6

[moonman]

Thanks King, can u explain 'Fs' and fs in this command

The book also tells to plot magnitude response by this

h770=[];
L=50;
fs=12000;
fb=770;
for n=1:L
    h770(n)=(2/L)*cos(2*pi*fb*n/fs);
end

fs=8000;
ww=0:(pi/256):pi;
ff=ww/(2*pi)*fs;
H=freqz(h770,1,ww);
plot(ff,abs(H));
grid on;

Although the shape is same for both plots but magnitude varies

[Wayne King]

The big difference is that you have not plotted yours in dB

plot(ff,abs(H));

% if you plot

plot(ff,20*log10(abs(H)));

You'll see. I think it is much more common to plot these in dB.

Answer 7

[Wayne King]

fvtool() is doing that under the hood and much more, I question why your book constructs a frequency axis in angular frequencies, when in this application Hertz makes much more sense:

    [H,F] = freqz(h,1,[],fs);
    plot(F./1000,20*log10(abs(H)));
    grid on;
    xlabel('kHz'); ylabel('Magnitude-squared (dB)');

Again, I would recommend you avoid a for loop to calculate your FIR filter coefficients and just use the vector approach I showed above.

Answer 8

[moonman]

Yes u are absolutely right I question why your book constructs a frequency axis in angular frequencies, when in this application Hertz makes much more sense:

My last query , the question says to

indicate the location of each of the DTMF frequencies(697,770,852,941,1209,1336 and 1477 Hz) on the plot generated by

h770=[];
L=50;
fs=12000;
fb=770;
for n=1:L
    h770(n)=(2/L)*cos(2*pi*fb*n/fs);
end
fs=8000;
ww=0:(pi/256):pi;
ff=ww/(2*pi)*fs;
H=freqz(h770,1,ww);
plot(ff,abs(H));
grid on;

How can i do that, kindly make me clear

Answer 9

[Wayne King]

    [H,F] = freqz(h,1,[],fs);
    plot(F,20*log10(abs(H)));
    grid on;
    xlabel('Hz'); ylabel('Magnitude-squared (dB)');
    set(gca,'xtick',[697,770,852,941,1209,1336,1477]); 
    set(gca,'xlim',[500 2000]);

Note that you really need to limit your x-axis in order to make the labels reasonable. If you use the whole frequency axis from 0 to the Nyquist, they bunch up.

Answer 10

[moonman]

Wonderful Help Wayne King Thanks a lot for making me understand this issue

Answer 11

[moonman]

Wayne King u wrote this code for me

   [H,F] = freqz(h,1,[],fs);
      plot(F,20*log10(abs(H)));
      grid on;
      xlabel('Hz'); ylabel('Magnitude-squared (dB)');
      set(gca,'xtick',[697,770,852,941,1209,1336,1477]); 
      set(gca,'xlim',[500 2000]);

Can u kindly explain me why u wrote Magnitude-Squared (db) in y axis, why squared

[moonman]

sir kindly explain why u used Magnitude Squared in y-axis

[Wayne King]

It is often desirable to obtain power estimates, which is the magnitude squared.

[moonman]

Can u kindly refer some webpage which further explain this concept.

Thanks a lot

[praveen]

kindly send total program in a single page of dtmf decoder using fir filters try to send quickleyyyyyy please

[Jan]

@praveen: This request is such unfriendly, that it is funny, that you hope to be successful. It is definitely your turn to pick from these many pieces of code a program, that is working for you. Trying to push us to do your work quickly is really the wrong approach.