Matlab Summation equation for calculating impulse response of FBMC

Hello guys,

I am trying to solve this equation to calculate the impulse response of a filter bank carrier but the graph is not exactly same to what i should get.

Kindly help!

M = 64; %Number of subcarriers

K = 4; %overlapping factor

y=[1 0.97195983 sqrt(2)/2 0.23514695 ]

T= M;

t = linspace (0,100,100)

for k = 1:3

result =(y(1,1)^2* cos(2*pi*k*t/(K*T)) + y(1,2)^2*cos(2*pi*k*t/(K*T))+ y(1,3)^2* cos(2*pi*k*t/(K*T)+ y(1,4)^2* cos(2*pi*k*t/(K*T))));

k

end

result = 1+2*result

plot (result)

Answer 1

[AJAY CHANDRA DORAGARI]

i understand that you are trying to add harmonics (fundamental upto third ) by computing harmonics at some instants of time and trying to get plot of it

M = 64; %Number of subcarriers

K = 4; %overlapping factor

y=[1 0.97195983 sqrt(2)/2 0.23514695 ]

T= M;

l=K*T; % simplify the equation instead of typing it every time

for k=1:3

for t=0:100

h(k,t+1)=((y(1)^2)*cos((2*k*pi*t)/l))+((y(2)^2)*cos((2*k*pi*t)/l))+((y(3)^2)*cos((2*k*pi*t)/l))+((y(4)^2)*cos((2*k*pi*t)/l))

%here im using t+1 because t=0 cannot be used since indexing must be a positve non zero interger i could use logic(t) but im too lazy

end

end

a=1+(2*sum(h))

plot(a)

xlabel('h(t)')

ylabel('t')

is this the graph you are expecting?

even code can be improved

[AJAY CHANDRA DORAGARI]

i would even like to know whether you are trying to get discrete or a continous one

i assumed that to be continous time function

since you have specified it as a function of 't' rather than 'n'

[Sara Nasir]

Hello Ajay,

I am trying to plot the prototype filter in time and frequency domain for which I have to solve the respective formulas. I am just starting of with Matlab coding since it is the part of my project. The graph I should get is something like this.

[AJAY CHANDRA DORAGARI]

I could try if you could attach a file about prototype filter and the mathematical relations about response you want sorry im not an electronics and communications student iam an electrical student but I had a subject signals and systems in my course,so I can try it.I tried searching about prototype filter but couldnt understand the terms associated with the formula.If you could explain the terms I can help you to write the code

[Sara Nasir]

okay Ajay I have to plot the frequency response and the impulse response which I am trying to get but the graphs are a little bit different from what is expected.

I have to plot h(t) and H(f) in the chapter 3 of the paper.

[Sara Nasir]

Meanwhile I will do the changes in the code as you have mentioned above.

[Sara Nasir]

Chapter 3 is the topic. This filter known as PHYDYAS filter. I have simulate it in time and frequency domain (in Matlab) using the formulas given for K=overlapping factor=4, M= number of subcarrier = 64 to reproduce figure 6 and Figure 7 . For this I have to use formula of h2(t) and H2(f) in page number 10.

[AJAY CHANDRA DORAGARI]

so did you able to get it?

[Sara Nasir]

Its different.

[AJAY CHANDRA DORAGARI]

no i made a bluder in code so im trying to rectify it

[Sara Nasir]

ok

[AJAY CHANDRA DORAGARI]

i have a doubt consider the h(0) instant calculate the value at h(0)

substituting t=0 and k=1,2,3 as per the given formula

h(0)=1+2*(H1*cos(2*pi*0)+H2*cos(4*pi*0)+H3*cos(6*pi*0))

this would be h(0)=1+2(h1+h2+h3)=approx 7

nd in the time response it is given as amplitude =0 at t=0 instant

[Sara Nasir]

yes in graph at t=0, Amp is 0

[AJAY CHANDRA DORAGARI]

h(0) calculated manually and h(0) in the graph must match if not the formula or the graph or our interpretation of formula must be wrong

tell me Hk in formula are the values given in the table 1 of page 7 ?

and the value of T is M?

[Sara Nasir]

yes the values are from table 1 on page 7.

K=4 and M=64

T is the duration of multicarrier symbol and it is M.

I can also share two forums on which code is given but the formulas are different there.

https://fr.mathworks.com/matlabcentral/fileexchange/42271-fbmc-ofdm-prototype-filter-comparison

https://fr.mathworks.com/matlabcentral/fileexchange/45100-fbmc-modulation-demodulation

[AJAY CHANDRA DORAGARI]

got it change the formula of time response

it is h(t)=1+summation[((-1)^k)*2*Hk*cos(2*pi*kt/KT)

[AJAY CHANDRA DORAGARI]

there they missed the term (-1)^k ill update the code give me a min

[Sara Nasir]

okay from where did u find the formula?

[AJAY CHANDRA DORAGARI]

correct it the formula which they used in the first document which you sent

M = 256; %Number of subcarriers

K = 4; %overlapping factor

y=[1 0.97195983 sqrt(2)/2 0.23514695 ]

T= M;

l=K*T; % simplify the equation instead of typing it every time

t=0:1000;

h=1+2*((-y(1)*cos((2*pi.*t)/l))+(y(2)*cos((4*pi.*t)/l))+(-y(3)*cos((6*pi.*t)/l)));

plot(h)

ylabel('h(t)')

xlabel('t')

grid on

[AJAY CHANDRA DORAGARI]

in this code i didnt use the for loop instead used element wise operation

use element wise operation wherever it is possible it is faster than the loop operations

[Sara Nasir]

okay it seems easier now

[AJAY CHANDRA DORAGARI]

formula i have seen it first at the ieee paper but i dont know which paper it was then alternatingly added and substracted the terms Hk after trial and error i got it and then cross checked it with the graph and it made sense. Some terms must be substracted to get h(0) value to zero then i remembered generally the terms in fft analysis are alternating with positve and negative signs

[Sara Nasir]

This seems to be the fromula. Thank you

[Sara Nasir]

but you know for M= 64 it is giving the same result

[AJAY CHANDRA DORAGARI]

i didnt get you can send me the plot

[Sara Nasir]

M = 64; %Number of subcarriers

K = 4; %overlapping factor

y=[1 0.97195983 sqrt(2)/2 0.23514695 ]

T= M;

l=K*T; % simplify the equation instead of typing it every time

t=0:1000;

h=1+2*((-y(1)*cos((2*pi.*t)/l))+(y(2)*cos((4*pi.*t)/l))+(-y(3)*cos((6*pi.*t)/l)));

plot(h)

ylabel('h(t)')

xlabel('t')

grid on

[Sara Nasir]

I think we should to 256 carriers

[AJAY CHANDRA DORAGARI]

no !

for 64 time period take upto t= 200

if you take more than 200 it will repeat again it will appear to be some sine wave where you cant make any distinction

[AJAY CHANDRA DORAGARI]

for 256

you take upto t= 1000 points

just play with it change the values and see the results check for different values of M

[Sara Nasir]

No the vaules are different. So it was the formula that was wrong.

[AJAY CHANDRA DORAGARI]

im sorry my bad not time period in the previous post

take the value of t=200 for M=64

t=1000 for M=256 (it is small t not capital)

[Sara Nasir]

got it! Thankyou

And in the similar way i should calculate the freq response by the formula H(f). No need of for loop?

[AJAY CHANDRA DORAGARI]

yes it was the formula which was wrong and ofcourse the code which i used in the first post was wrong too which was a blunder i didnt cross check the formula which you had mentioned .i have taken the expression from your code which you have posted i made a mistake and realised it but in the end i could rectify it

[Sara Nasir]

yes the formula is the paper was wrong

[Sara Nasir]

i will try to solve H(f) now and if I have any problem I will ask for help

[AJAY CHANDRA DORAGARI]

you could getaway maybe this time but even without for loop for frequency analysis the expression might be bit clumpsy because it contains six terms which can be simplified but im no expert i guess you could use it but sometimes you want have any other choice rather than using for loop

[AJAY CHANDRA DORAGARI]

i should thanq sara for keeping me busy during this pandemic time i could learn something about phydyas filter and rewind my coding even i started learning matlab just 2 months back for my masters project please do me a favour dont delete this question so i could revisit someday it might help for my own filter problem

[Sara Nasir]

I am trying to do with for loop because it is difficult to keep track of six values.

Well I am getting error this time during implementing it.

[AJAY CHANDRA DORAGARI]

can you tell me the error

[Sara Nasir]

Can you go through the code please!

There is no error now but nothing on graph.

M = 64; %Number of subcarriers

K = 4; %overlapping factor

H0=1;

H1=0.97195983;

H2=sqrt(2)/2;

H3=0.23514695;

H4=-0.97195983;

H5=-sqrt(2)/2;

H6=-0.23514695;

a=M*K;

f=-1/M:0.001:1/M;

for k=-3:3

for f=-1/M:0.001:1/M

b = f-k/a;

Num = sin(pi*(b)*a);

Den = a* sin(pi*(b));

product=Num./Den;

H = H6*product + H5*product + H4*product + H3*product + H2*product + H1*product + H0*product;

end

end

ab=1+(2*sum(H))

figure(1)

plot (ab)

xlabel('f')

ylabel('H(f)')

[AJAY CHANDRA DORAGARI]

i have a doubt

H-1 and H1 are both same the terms used in the frequency response formula

ill try to write my own code first then ill see your code