24 views (last 30 days)

Show older comments

HI i have a matlab code but i am unable to find the mistake. i have attched the code.

Basically what i am unable to achieve is based on the Quality factor formula Q=f/f2-f1 we supposed to have the same quality factor in the input values alpha = 2*pi/Q. so, When, I change Amplitude or the frequency the Q factor goes far away from the right Q factor. In, unfortunately, I think I got some errors in the dB scale or something else. I must have Q factor in the output as the same Q factor in the input. If wouldn't mind I submit the Matlab code take look at and fix the problem. I cant figure out where is the error in the code that is messing up with this.

- Input:
- Digital Signal:
- F= 100 Hz
- Q= 500
- alpha= (f*pi)/Q
- A= 10

Result:

Q= fmax/f2-f1

Q= 237.5527

- Digital Signal:
- F= 100 Hz
- Q= 500
- A= 100

Result:

Q= 99.4541

Any help is appreciated. Thanks in advance.

clear all

clc

close all

fs= 10000

t = (0:1/fs:(1-(1/fs)));

f= 4000

alpha= (f*pi)/500

y = 10*exp(-alpha*t).*sin(2*pi*f*t);

figure

plot(t,y)

N1=10000

S =fft(y,N1);

F1=S(1:N1/2); %half of spectra

PF1=2*F1.*conj(F1)/(fs*N1); %Power spectrum density per 1 kHz

E=sum(PF1)%energy of signal V^2*T

Xdb = 20*log10(S);

LPF1=10*log10(PF1); %Power spectrum in dB scale

w=(1:N1/2);

LP=LPF1(1:N1/2);

w1=fs*w/N1;

figure

plot(w1,LP)

BB= LP

CC= w1

E1 = max(LP)

t1= max(t)

talph1= E / -t1

%alpha1= max(BB*20*log(e))

[max_Z, max_index]=max(BB) %maximum value of impedance

threedb=max_Z*sqrt(2)/2; %the 3db point

[db_Z,index_db] = min(abs(BB.'-threedb)-max_Z)

%db3idx = zci(pwr-hpp);

[pk,loc] = findpeaks(BB,'Npeaks',1,'SortStr','descend');

%db3c= threedb-max_Z

ofst= 12

for k1 = 1:length(pk)

varmtx=[(BB(loc(k1)-ofst:loc(k1)));CC(loc(k1)-ofst:loc(k1));(BB(loc(k1):loc(k1)+ofst));CC(loc(k1):loc(k1)+ofst)]

dBpts(k1,1:2) = interp1(varmtx(1,:), varmtx(2,:), -(db_Z), 'linear','extrap');

dBpts(k1,3:4) = interp1(varmtx(3,:),varmtx(4,:), -(db_Z), 'linear','extrap');

end

figure(3)

plot(CC, BB)

hold on

for k1 = 1:length(pk)

plot(dBpts(k1,:), -(db_Z), 'r+')

end

hold off

grid

F_1 = dBpts(k1,1:2)

F_1 = F_1(1)

F_2 = dBpts(k1,3:4)

F_2 =F_2(1)

BW = (F_2)-(F_1)

Q = abs((loc-1)*(1/BW))

David Goodmanson
on 26 Feb 2021

Hello Abdullah

be aware that what you are calling alpha is often called alpha/2, but I don't think that had anything to do with the problem you are having. I stayed with your definition of alpha.

Everything appears to be fine up until the width calculation. I don't know whether you intended that w denote circular frequency, but in this case it is (regular) frequency. To avoid ambiguity I went with frequency and did not use circular frequency. The code below finds the widths and the resulting Q.

clc

close all

fs= 10000;

t = (0:1/fs:(1-(1/fs)));

f= 4000;

alpha= (f*pi)/500;

y = 10*exp(-alpha*t).*sin(2*pi*f*t);

figure(1)

plot(t,y); grid on

Q0 = pi*f/alpha

N1=10000;

S =fft(y,N1);

F1=S(1:N1/2); %half of spectra

PF1=2*F1.*conj(F1)/(fs*N1); %Power spectrum density per 1 kHz

% E=sum(PF1);%energy of signal V^2*T

% Xdb = 20*log10(S);

LPF1=10*log10(PF1); %Power spectrum in dB scale

LP=LPF1(1:N1/2);

fvec=fs*(1:N1/2)/N1; % frequency array

figure(2)

plot(fvec,LP); grid on

xlim([3990 4010])

% find half power points by interpolation of LP

[LPmax ind0] = max(LP);

LPhalf = LPmax - 10*log10(2); % half power points

indx = find(LP>LPhalf);

ind1 = [indx(1)-5:ind0]; % use a few more index points, 5 is somewhat arbitrary

f1 = interp1(LP(ind1),fvec(ind1),LPhalf,'spline');

ind2 = [ind0:indx(end)+5];

f2 = interp1(LP(ind2),fvec(ind2),LPhalf,'spline');

Q = f/(f2-f1)

Q0 = 500.0000

Q = 499.9917

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!