Thread Subject: Problem with "for"

I do not konw why but I have a problem with one of my "for" in my code.
Please help me!!! LOL

This is my code and my problem is in the 3rd "for" with "k=(N/2+2):N"
because the caculation of wc(k) give twice the number of wc I want. In fact
I obtain N values of wc while I need N/2 values and between 1 and N/2 all
the values are zero. I do not understand why. Yesterday when I did few
changes, I obtaines 3 times N values.
Thank you for your help

%%% vérifier longueur point%%%
n=1;
while 2^n<Zu.npts;
    n=n+1;
end;
corr=2^n;
plus=corr-Zu.npts;
ajout=zeros(plus,1);
new=vertcat(Zu.trace,ajout);
N=length(new);

% Transformée de fourier discrète du signal : a(t)-> A(w):
A=fft(new);
% Aa=A(1:N/2);
% Ab=A(N/2+1:N);

%Damping de 5% :
epsi=0.05;

% %Complex :
% i=sqrt(-1);

%%% Matrice vide %%%
Sd=[];

%%%pas en fréquence %%%%
pas=1/(N*Zu.delta);

%%% fréquence propre en Hz %%%
delta=10^(1/30);

for m=1:91;
    wo=0.1*(delta)^m;

    for j=1:N/2+1;
        w(j)=2*pi*pas*(j-1);
        %Fonction de transfert H en fonction de w, wo et epsi :
        H(j)=1/(w(j)^2-(2*i*(wo)*(epsi)*w(j))-wo^2);
        Xa(j)=A(j,1)*H(j);
    end


    for k=(N/2+2):N;
        wc(k)=-2*pi*pas*(N+1-k);
        %Fonction de transfert H en fonction de w, wo et epsi :
        Hc(k)=1/(wc(k)^2-(2*i*(wo)*(epsi)*wc(k))-wo^2);
        Xb(k)=A(k,1)*Hc(k);
    end


    X=vertcat(Xa,Xb);

     x=ifft(X);
     ma=max(x);
     Sd=[Sd;[ma wo]];
     S=vertcat(Sd);
end
    semilogx(S(:,2),S(:,1));
    xlabel ('frequency (Hz)','FontSize',10); ylabel
('displacement','FontSize',10);
    title('Displacement response spectra','FontSize',14);

"Estelle" <tahoser45@hotmail.com> wrote in message
<477b5bff$0$21142$7a628cd7@news.club-internet.fr>...
> I do not konw why but I have a problem with one of my "for"
in my code.
> Please help me!!! LOL
>
> This is my code and my problem is in the 3rd "for" with
"k=(N/2+2):N"
> because the caculation of wc(k) give twice the number of wc I
want. In fact
> I obtain N values of wc while I need N/2 values and between
1 and N/2 all
> the values are zero. I do not understand why. Yesterday when
I did few
> changes, I obtaines 3 times N values.
> Thank you for your help
>
> %%% vérifier longueur point%%%
> n=1;
> while 2^n<Zu.npts;
> n=n+1;
> end;
> corr=2^n;
> plus=corr-Zu.npts;
> ajout=zeros(plus,1);
> new=vertcat(Zu.trace,ajout);
> N=length(new);
>
> % Transformée de fourier discrète du signal : a(t)-> A(w):
> A=fft(new);
> % Aa=A(1:N/2);
> % Ab=A(N/2+1:N);
>
> %Damping de 5% :
> epsi=0.05;
>
> % %Complex :
> % i=sqrt(-1);
>
> %%% Matrice vide %%%
> Sd=[];
>
> %%%pas en fréquence %%%%
> pas=1/(N*Zu.delta);
>
> %%% fréquence propre en Hz %%%
> delta=10^(1/30);
>
> for m=1:91;
> wo=0.1*(delta)^m;
>
> for j=1:N/2+1;
> w(j)=2*pi*pas*(j-1);
> %Fonction de transfert H en fonction de w, wo et epsi
 :
> H(j)=1/(w(j)^2-(2*i*(wo)*(epsi)*w(j))-wo^2);
> Xa(j)=A(j,1)*H(j);
> end
>
>
> for k=(N/2+2):N;
> wc(k)=-2*pi*pas*(N+1-k);
> %Fonction de transfert H en fonction de w, wo et epsi
 :
> Hc(k)=1/(wc(k)^2-(2*i*(wo)*(epsi)*wc(k))-wo^2);
> Xb(k)=A(k,1)*Hc(k);
> end
>
>
> X=vertcat(Xa,Xb);
>
> x=ifft(X);
> ma=max(x);
> Sd=[Sd;[ma wo]];
> S=vertcat(Sd);
> end
> semilogx(S(:,2),S(:,1));
> xlabel ('frequency (Hz)','FontSize',10); ylabel
> ('displacement','FontSize',10);
> title('Displacement response spectra','FontSize',14);
>
>
>
>
>

This example could illustrate your error:

a(3) = 1

The result will be

a = [0 0 1]

Because k gets values in the range [N/2+2),N], wc becomes a
vector with N elements.

Moreover, the soure code needs to be vectorized and avoid using
loop. Here is just an example which could replace the 2nd and
3rd loops

    w = 2*pi*pas*[0:N/2 -N/2+1:-1];
    H = 1./(w.^2-2*i*wo*epsi*w-wo^2);
    X = A.* H';

Anh Huy Phan
RIKEN - BSI

Hi,
Thanks for your help.
I just do not understand why you said that wc has N values. Because the
values are included into N/2 and N. I know that it could be evident for you
but not for me and I really want to understand the problem.
I try your code as well but it does not work.

I am sorry if I am a bit insistent but I really want to understand all
steps!

Thank you very much
Best regards
Estelle


"Huy " <phananhhuy@mathworks.com> a écrit dans le message de news:
flfr61$mtu$1@fred.mathworks.com...
> "Estelle" <tahoser45@hotmail.com> wrote in message
> <477b5bff$0$21142$7a628cd7@news.club-internet.fr>...
>> I do not konw why but I have a problem with one of my "for"
> in my code.
>> Please help me!!! LOL
>>
>> This is my code and my problem is in the 3rd "for" with
> "k=(N/2+2):N"
>> because the caculation of wc(k) give twice the number of wc I
> want. In fact
>> I obtain N values of wc while I need N/2 values and between
> 1 and N/2 all
>> the values are zero. I do not understand why. Yesterday when
> I did few
>> changes, I obtaines 3 times N values.
>> Thank you for your help
>>
>> %%% vérifier longueur point%%%
>> n=1;
>> while 2^n<Zu.npts;
>> n=n+1;
>> end;
>> corr=2^n;
>> plus=corr-Zu.npts;
>> ajout=zeros(plus,1);
>> new=vertcat(Zu.trace,ajout);
>> N=length(new);
>>
>> % Transformée de fourier discrète du signal : a(t)-> A(w):
>> A=fft(new);
>> % Aa=A(1:N/2);
>> % Ab=A(N/2+1:N);
>>
>> %Damping de 5% :
>> epsi=0.05;
>>
>> % %Complex :
>> % i=sqrt(-1);
>>
>> %%% Matrice vide %%%
>> Sd=[];
>>
>> %%%pas en fréquence %%%%
>> pas=1/(N*Zu.delta);
>>
>> %%% fréquence propre en Hz %%%
>> delta=10^(1/30);
>>
>> for m=1:91;
>> wo=0.1*(delta)^m;
>>
>> for j=1:N/2+1;
>> w(j)=2*pi*pas*(j-1);
>> %Fonction de transfert H en fonction de w, wo et epsi
> :
>> H(j)=1/(w(j)^2-(2*i*(wo)*(epsi)*w(j))-wo^2);
>> Xa(j)=A(j,1)*H(j);
>> end
>>
>>
>> for k=(N/2+2):N;
>> wc(k)=-2*pi*pas*(N+1-k);
>> %Fonction de transfert H en fonction de w, wo et epsi
> :
>> Hc(k)=1/(wc(k)^2-(2*i*(wo)*(epsi)*wc(k))-wo^2);
>> Xb(k)=A(k,1)*Hc(k);
>> end
>>
>>
>> X=vertcat(Xa,Xb);
>>
>> x=ifft(X);
>> ma=max(x);
>> Sd=[Sd;[ma wo]];
>> S=vertcat(Sd);
>> end
>> semilogx(S(:,2),S(:,1));
>> xlabel ('frequency (Hz)','FontSize',10); ylabel
>> ('displacement','FontSize',10);
>> title('Displacement response spectra','FontSize',14);
>>
>>
>>
>>
>>
>
> This example could illustrate your error:
>
> a(3) = 1
>
> The result will be
>
> a = [0 0 1]
>
> Because k gets values in the range [N/2+2),N], wc becomes a
> vector with N elements.
>
> Moreover, the soure code needs to be vectorized and avoid using
> loop. Here is just an example which could replace the 2nd and
> 3rd loops
>
> w = 2*pi*pas*[0:N/2 -N/2+1:-1];
> H = 1./(w.^2-2*i*wo*epsi*w-wo^2);
> X = A.* H';
>
> Anh Huy Phan
> RIKEN - BSI
>
>
>
>
>
>
>
>
>

Assume that N = 8, you type this line in the command window:

  wc(N/2+2) = 1
  
-> wc =

     0 0 0 0 0 1

You get the vector wc with 8/2+2 = 6 elements

And then, type this line

  wc(N) = 1

-> wc =
  
     0 0 0 0 0 1 0 1

wc has 8 elements.

And the modified source code is as below

...
for m=1:91;
    wo=0.1*(delta)^m;
    
    %%%
    w = 2*pi*pas*[0:N/2 -N/2+1:-1];
    H = 1./(w.^2-2*i*wo*epsi*w-wo^2);
    X = A.* H';
    %%%

    x=ifft(X);
    ma=max(x);
    Sd=[Sd;[ma wo]];
    S=vertcat(Sd);
end
...

By the ways, there is another way without using for loop.

Anh Huy Phan
RIKEN - BSI


"Estelle" <tahoser45@hotmail.com> wrote in message
<477b8cdd$0$21146$7a628cd7@news.club-internet.fr>...
> Hi,
> Thanks for your help.
> I just do not understand why you said that wc has N values.
Because the
> values are included into N/2 and N. I know that it could be
evident for you
> but not for me and I really want to understand the problem.
> I try your code as well but it does not work.
>
> I am sorry if I am a bit insistent but I really want to
understand all
> steps!
>
> Thank you very much
> Best regards
> Estelle
>
>
> "Huy " <phananhhuy@mathworks.com> a écrit dans le message de
news:
> flfr61$mtu$1@fred.mathworks.com...
> > "Estelle" <tahoser45@hotmail.com> wrote in message
> > <477b5bff$0$21142$7a628cd7@news.club-internet.fr>...
> >> I do not konw why but I have a problem with one of my "for"
> > in my code.
> >> Please help me!!! LOL
> >>
> >> This is my code and my problem is in the 3rd "for" with
> > "k=(N/2+2):N"
> >> because the caculation of wc(k) give twice the number of
wc I
> > want. In fact
> >> I obtain N values of wc while I need N/2 values and between
> > 1 and N/2 all
> >> the values are zero. I do not understand why. Yesterday
when
> > I did few
> >> changes, I obtaines 3 times N values.
> >> Thank you for your help
> >>
> >> %%% vérifier longueur point%%%
> >> n=1;
> >> while 2^n<Zu.npts;
> >> n=n+1;
> >> end;
> >> corr=2^n;
> >> plus=corr-Zu.npts;
> >> ajout=zeros(plus,1);
> >> new=vertcat(Zu.trace,ajout);
> >> N=length(new);
> >>
> >> % Transformée de fourier discrète du signal : a(t)-> A(w):
> >> A=fft(new);
> >> % Aa=A(1:N/2);
> >> % Ab=A(N/2+1:N);
> >>
> >> %Damping de 5% :
> >> epsi=0.05;
> >>
> >> % %Complex :
> >> % i=sqrt(-1);
> >>
> >> %%% Matrice vide %%%
> >> Sd=[];
> >>
> >> %%%pas en fréquence %%%%
> >> pas=1/(N*Zu.delta);
> >>
> >> %%% fréquence propre en Hz %%%
> >> delta=10^(1/30);
> >>
> >> for m=1:91;
> >> wo=0.1*(delta)^m;
> >>
> >> for j=1:N/2+1;
> >> w(j)=2*pi*pas*(j-1);
> >> %Fonction de transfert H en fonction de w, wo et
epsi
> > :
> >> H(j)=1/(w(j)^2-(2*i*(wo)*(epsi)*w(j))-wo^2);
> >> Xa(j)=A(j,1)*H(j);
> >> end
> >>
> >>
> >> for k=(N/2+2):N;
> >> wc(k)=-2*pi*pas*(N+1-k);
> >> %Fonction de transfert H en fonction de w, wo et
epsi
> > :
> >> Hc(k)=1/(wc(k)^2-(2*i*(wo)*(epsi)*wc(k))-wo^2);
> >> Xb(k)=A(k,1)*Hc(k);
> >> end
> >>
> >>
> >> X=vertcat(Xa,Xb);
> >>
> >> x=ifft(X);
> >> ma=max(x);
> >> Sd=[Sd;[ma wo]];
> >> S=vertcat(Sd);
> >> end
> >> semilogx(S(:,2),S(:,1));
> >> xlabel ('frequency (Hz)','FontSize',10); ylabel
> >> ('displacement','FontSize',10);
> >> title('Displacement response spectra','FontSize',14);
> >>
> >>
> >>
> >>
> >>
> >
> > This example could illustrate your error:
> >
> > a(3) = 1
> >
> > The result will be
> >
> > a = [0 0 1]
> >
> > Because k gets values in the range [N/2+2),N], wc becomes
a
> > vector with N elements.
> >
> > Moreover, the soure code needs to be vectorized and avoid
using
> > loop. Here is just an example which could replace the 2nd
and
> > 3rd loops
> >
> > w = 2*pi*pas*[0:N/2 -N/2+1:-1];
> > H = 1./(w.^2-2*i*wo*epsi*w-wo^2);
> > X = A.* H';
> >
> > Anh Huy Phan
> > RIKEN - BSI
> >
> >
> >
> >
> >
> >
> >
> >
> >
>
>

OK thank you!
"Huy " <phananhhuy@mathworks.com> a écrit dans le message de news:
flg5l8$917$1@fred.mathworks.com...
> Assume that N = 8, you type this line in the command window:
>
> wc(N/2+2) = 1
>
> -> wc =
>
> 0 0 0 0 0 1
>
> You get the vector wc with 8/2+2 = 6 elements
>
> And then, type this line
>
> wc(N) = 1
>
> -> wc =
>
> 0 0 0 0 0 1 0 1
>
> wc has 8 elements.
>
> And the modified source code is as below
>
> ...
> for m=1:91;
> wo=0.1*(delta)^m;
>
> %%%
> w = 2*pi*pas*[0:N/2 -N/2+1:-1];
> H = 1./(w.^2-2*i*wo*epsi*w-wo^2);
> X = A.* H';
> %%%
>
> x=ifft(X);
> ma=max(x);
> Sd=[Sd;[ma wo]];
> S=vertcat(Sd);
> end
> ...
>
> By the ways, there is another way without using for loop.
>
> Anh Huy Phan
> RIKEN - BSI
>
>
> "Estelle" <tahoser45@hotmail.com> wrote in message
> <477b8cdd$0$21146$7a628cd7@news.club-internet.fr>...
>> Hi,
>> Thanks for your help.
>> I just do not understand why you said that wc has N values.
> Because the
>> values are included into N/2 and N. I know that it could be
> evident for you
>> but not for me and I really want to understand the problem.
>> I try your code as well but it does not work.
>>
>> I am sorry if I am a bit insistent but I really want to
> understand all
>> steps!
>>
>> Thank you very much
>> Best regards
>> Estelle
>>
>>
>> "Huy " <phananhhuy@mathworks.com> a écrit dans le message de
> news:
>> flfr61$mtu$1@fred.mathworks.com...
>> > "Estelle" <tahoser45@hotmail.com> wrote in message
>> > <477b5bff$0$21142$7a628cd7@news.club-internet.fr>...
>> >> I do not konw why but I have a problem with one of my "for"
>> > in my code.
>> >> Please help me!!! LOL
>> >>
>> >> This is my code and my problem is in the 3rd "for" with
>> > "k=(N/2+2):N"
>> >> because the caculation of wc(k) give twice the number of
> wc I
>> > want. In fact
>> >> I obtain N values of wc while I need N/2 values and between
>> > 1 and N/2 all
>> >> the values are zero. I do not understand why. Yesterday
> when
>> > I did few
>> >> changes, I obtaines 3 times N values.
>> >> Thank you for your help
>> >>
>> >> %%% vérifier longueur point%%%
>> >> n=1;
>> >> while 2^n<Zu.npts;
>> >> n=n+1;
>> >> end;
>> >> corr=2^n;
>> >> plus=corr-Zu.npts;
>> >> ajout=zeros(plus,1);
>> >> new=vertcat(Zu.trace,ajout);
>> >> N=length(new);
>> >>
>> >> % Transformée de fourier discrète du signal : a(t)-> A(w):
>> >> A=fft(new);
>> >> % Aa=A(1:N/2);
>> >> % Ab=A(N/2+1:N);
>> >>
>> >> %Damping de 5% :
>> >> epsi=0.05;
>> >>
>> >> % %Complex :
>> >> % i=sqrt(-1);
>> >>
>> >> %%% Matrice vide %%%
>> >> Sd=[];
>> >>
>> >> %%%pas en fréquence %%%%
>> >> pas=1/(N*Zu.delta);
>> >>
>> >> %%% fréquence propre en Hz %%%
>> >> delta=10^(1/30);
>> >>
>> >> for m=1:91;
>> >> wo=0.1*(delta)^m;
>> >>
>> >> for j=1:N/2+1;
>> >> w(j)=2*pi*pas*(j-1);
>> >> %Fonction de transfert H en fonction de w, wo et
> epsi
>> > :
>> >> H(j)=1/(w(j)^2-(2*i*(wo)*(epsi)*w(j))-wo^2);
>> >> Xa(j)=A(j,1)*H(j);
>> >> end
>> >>
>> >>
>> >> for k=(N/2+2):N;
>> >> wc(k)=-2*pi*pas*(N+1-k);
>> >> %Fonction de transfert H en fonction de w, wo et
> epsi
>> > :
>> >> Hc(k)=1/(wc(k)^2-(2*i*(wo)*(epsi)*wc(k))-wo^2);
>> >> Xb(k)=A(k,1)*Hc(k);
>> >> end
>> >>
>> >>
>> >> X=vertcat(Xa,Xb);
>> >>
>> >> x=ifft(X);
>> >> ma=max(x);
>> >> Sd=[Sd;[ma wo]];
>> >> S=vertcat(Sd);
>> >> end
>> >> semilogx(S(:,2),S(:,1));
>> >> xlabel ('frequency (Hz)','FontSize',10); ylabel
>> >> ('displacement','FontSize',10);
>> >> title('Displacement response spectra','FontSize',14);
>> >>
>> >>
>> >>
>> >>
>> >>
>> >
>> > This example could illustrate your error:
>> >
>> > a(3) = 1
>> >
>> > The result will be
>> >
>> > a = [0 0 1]
>> >
>> > Because k gets values in the range [N/2+2),N], wc becomes
> a
>> > vector with N elements.
>> >
>> > Moreover, the soure code needs to be vectorized and avoid
> using
>> > loop. Here is just an example which could replace the 2nd
> and
>> > 3rd loops
>> >
>> > w = 2*pi*pas*[0:N/2 -N/2+1:-1];
>> > H = 1./(w.^2-2*i*wo*epsi*w-wo^2);
>> > X = A.* H';
>> >
>> > Anh Huy Phan
>> > RIKEN - BSI
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>> >
>>
>>
>