Question Fourier Method, DFT and irregular sampling
in [Matlab]
|
Prev: Symbol Line Width in Legend
Next: ode solver options
From: kk KKsingh on 31 May 2010 14:57 I read one of the paper ! It was about reconstrution of a simple sin signal ! Greag please do put your comments on it ! Objective : To reconstruct the simple sin wave through Dujhdham and Schonvillie algorithm Methodology : 1. Take DFT 2. Apply weights 3. Multiply with the inverse kernal 4. Simple ifft will reconstruct the samples back My code ( Some thing is wrong) clear all; close all; N=200; dt=.002; t=0:dt:dt*(N-1); signal=sin(2*pi*10*t); %________--------------------------------------------------- %_choosing percentage of decimation__ % Irregular interpolation % Randomizing the offset locations kn1=randperm(N); % Randomly eliminating some traces percent=0.25 kn2=kn1(1:ceil(percent*N)); % sorting random chosen traces in ascending order z=[sort(kn2)]; ------------------------------------------------------------------------- y1=signal; t1=t; % Finding the indices of known components y1(z)=[]; t1(z)=[] % Making the weight vector ! Period=max(t1)-min(t1); deltax_n=zeros(N2,1); deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; for k=2:N2-1; deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); end %NDFT with weights N2=length(t1); X2 = zeros(size(y1)); for k=1:N2 freq = (-1j*2*pi*(k-1))/N2; Wnk = exp(freq.*(0:length(t1)-1)); X2(k)=sum(y1.*Wnk*deltax_n(k)); end % Making the weight vector ! Period=max(t1)-min(t1); deltax_n=zeros(N2,1); deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; for k=2:N2-1; deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); end %DFT with weights N2=length(t1); X2 = zeros(size(y1)); M=N2/2; %(Author told me to always take more N and less M to make the reconstruction possible) for k=1:M freq = (-1j*2*pi*(k-1))/N2; Wnk = exp(freq.*(0:length(t1)-1)); X2(k)=sum(y1.*Wnk*deltax_n(k)); end %Now constructing the inverse kernal %It should look like this finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk) * X2; %Although expression is like finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk+KI) * X2; %but author stress that if i take more number of samples than fourier cofecients than i %dont need eye more over eye is just to suppress noise in any signal result=ifft(final spectrum); %I am doing some thing wrong here ! Author of the paper told me it should work but its not working for me ! Please guide me ! %If my question is not complete please let me know
From: kk KKsingh on 31 May 2010 15:12 > I read one of the paper ! It was about reconstrution of a simple sin signal ! Greag please do put your comments on it ! > > Objective : To reconstruct the simple sin wave through Dujhdham and Schonvillie algorithm > > Methodology : 1. Take DFT > 2. Apply weights > 3. Multiply with the inverse kernal > 4. Simple ifft will reconstruct the samples back > > > My code ( Some thing is wrong) %Aki Corrected code clear all; close all; N=200; dt=.002; t=0:dt:dt*(N-1); signal=sin(2*pi*10*t); %save unisin.dat signal -ascii %Plotting the sin signal figure(1) hold on stem(t,signal,'r') plot(t,signal) xlabel('time','fontsize',16) ylabel('Amplitude','fontsize',16) set(gca,'fontsize',16); title('Sin signal'); set(gcf,'color','w'); set(gca,'linewidth',3) %Making Sample Missing from the Simple Sin Signal %________--------------------------------------------------- %_choosing percentage of decimation__ % Irregular interpolation % Randomizing the offset locations kn1=randperm(N); % Randomly eliminating some traces percent=0.25 kn2=kn1(1:ceil(percent*N)); % sorting random chosen traces in ascending order z=[sort(kn2)]; %------------------------------------------------------------------------- y1=signal; t1=t; % %h1=[1:N]; [%a,b]=wcommon(z,h1); %H=find(b-1); y1(z)=[]; t1(z)=[]; %save nonunisin.dat y1 -ascii %fid=fopen('nonunisin.dat','w'); %printf(fid,'%d\n',y1); %fclose(fid); %save nonunitime.dat t1 -ascii %figure(2) %hold on %plot(t,signal) %stem(t(H),signal(H),'r') %hold off %xlabel('time','fontsize',16) %ylabel('Amplitude','fontsize',16) %set(gca,'fontsize',16); %itle('Decimated signal'); %set(gcf,'color','w'); %set(gca,'linewidth',3) %Writing code for DFT ( Discrete Fourier Transform) %There are no weights N1=length(t1); X1 = zeros(size(y1)); M=N1/2; for k=1:M freq = (-1j*2*pi*(k-1))/N1; Wnk = exp(freq.*(0:length(t1)-1)); X1(k)=sum(y1.*Wnk); end if mod(M,2)==0 k_dft=-N1/2:N1/2-1; % N even else k_dft=-(M-1)/2:(M-1)/2; % N odd end frqaxs_dft=1/dt/N*k_dft; figure(4) plot(frqaxs_dft,fftshift(abs(X1))) xlabel('frequency','fontsize',16) ylabel('Power','fontsize',16) set(gca,'fontsize',16); title('Energy of non uniform sin signal (Positive)'); set(gcf,'color','w'); %plotting the spectrum for DFT with Weights ! N2=length(y1) M2=N2/2; if mod(M2,2)==0 k_NDFT=-M2/2:M2/2-1; else k_NDFT=-(M2-1)/2:(M2-1)/2 end freq_NDFT=1/dt/M2*k_NDFT; % Making the weight vector ! Period=max(t1)-min(t1); deltax_n=zeros(N2,1); deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; deltax_n(N2,:)=((t1(1)+Period)-(t1(N2-1)))/2; for k=2:N2-1; deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); end %NDFT with weights N2=length(t1); X2 = zeros(size(y1)); M2=N2/2; for k=1:M2 freq = (-1j*2*pi*(k-1))/N2; Wnk = exp(freq.*(0:length(t1)-1)); X2(k)=sum(y1.*Wnk*deltax_n(k)); end %Contructiing a Fourier kernal %> %Now constructing the inverse kernal > %It should look like this finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk) * X2; > %Although expression is like finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk+KI) * X2; > %but author stress that if i take more number of samples than fourier coffecients than %dont need eye more it is just to suppress noise in any signal i:e for stabilization of %diagnonal > > > > result=ifft(final spectrum); > > "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu10pv$5dp$1(a)fred.mathworks.com>... > I read one of the paper ! It was about reconstrution of a simple sin signal ! Greag please do put your comments on it ! > > Objective : To reconstruct the simple sin wave through Dujhdham and Schonvillie algorithm > > Methodology : 1. Take DFT > 2. Apply weights > 3. Multiply with the inverse kernal > 4. Simple ifft will reconstruct the samples back > > > My code ( Some thing is wrong) > > clear all; > close all; > N=200; > dt=.002; > t=0:dt:dt*(N-1); > signal=sin(2*pi*10*t); > > > %________--------------------------------------------------- > %_choosing percentage of decimation__ > % Irregular interpolation > % Randomizing the offset locations > kn1=randperm(N); > % Randomly eliminating some traces > percent=0.25 > kn2=kn1(1:ceil(percent*N)); > % sorting random chosen traces in ascending order > z=[sort(kn2)]; > > > ------------------------------------------------------------------------- > y1=signal; t1=t; > > % Finding the indices of known components > > y1(z)=[]; > t1(z)=[] > % Making the weight vector ! > > Period=max(t1)-min(t1); > deltax_n=zeros(N2,1); > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > for k=2:N2-1; > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > end > > %NDFT with weights > > > N2=length(t1); > X2 = zeros(size(y1)); > > for k=1:N2 > freq = (-1j*2*pi*(k-1))/N2; > Wnk = exp(freq.*(0:length(t1)-1)); > X2(k)=sum(y1.*Wnk*deltax_n(k)); > end > > % Making the weight vector ! > > Period=max(t1)-min(t1); > deltax_n=zeros(N2,1); > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > for k=2:N2-1; > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > end > > %DFT with weights > > > N2=length(t1); > X2 = zeros(size(y1)); > M=N2/2; %(Author told me to always take more N and less M to make the reconstruction possible) > for k=1:M > freq = (-1j*2*pi*(k-1))/N2; > Wnk = exp(freq.*(0:length(t1)-1)); > X2(k)=sum(y1.*Wnk*deltax_n(k)); > end > > %Now constructing the inverse kernal > %It should look like this finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk) * X2; > %Although expression is like finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk+KI) * X2; > %but author stress that if i take more number of samples than fourier cofecients than i %dont need eye more over eye is just to suppress noise in any signal > > > > result=ifft(final spectrum); > > > > %I am doing some thing wrong here ! Author of the paper told me it should work but its not working for me ! Please guide me ! > > > %If my question is not complete please let me know
From: kk KKsingh on 31 May 2010 15:18 oops another mistake inverse kernal which we multiplying with DFT spcetrum is inv(Forward operater * W*Backward operater).... or inv(Forward operater * W*Backward operater+KI).... Thanks "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu11m3$15t$1(a)fred.mathworks.com>... > > I read one of the paper ! It was about reconstrution of a simple sin signal ! Greag please do put your comments on it ! > > > > Objective : To reconstruct the simple sin wave through Dujhdham and Schonvillie algorithm > > > > Methodology : 1. Take DFT > > 2. Apply weights > > 3. Multiply with the inverse kernal > > 4. Simple ifft will reconstruct the samples back > > > > > > My code ( Some thing is wrong) > > %Aki Corrected code > clear all; > close all; > N=200; > dt=.002; > t=0:dt:dt*(N-1); > signal=sin(2*pi*10*t); > > %save unisin.dat signal -ascii > > > %Plotting the sin signal > > figure(1) > hold on > stem(t,signal,'r') > plot(t,signal) > xlabel('time','fontsize',16) > ylabel('Amplitude','fontsize',16) > set(gca,'fontsize',16); > title('Sin signal'); > set(gcf,'color','w'); > set(gca,'linewidth',3) > > > %Making Sample Missing from the Simple Sin Signal > > %________--------------------------------------------------- > %_choosing percentage of decimation__ > % Irregular interpolation > % Randomizing the offset locations > kn1=randperm(N); > % Randomly eliminating some traces > percent=0.25 > kn2=kn1(1:ceil(percent*N)); > % sorting random chosen traces in ascending order > z=[sort(kn2)]; > > > > %------------------------------------------------------------------------- > y1=signal; t1=t; > > % > %h1=[1:N]; > [%a,b]=wcommon(z,h1); > %H=find(b-1); > y1(z)=[]; > t1(z)=[]; > > %save nonunisin.dat y1 -ascii > %fid=fopen('nonunisin.dat','w'); > %printf(fid,'%d\n',y1); > %fclose(fid); > %save nonunitime.dat t1 -ascii > > %figure(2) > %hold on > %plot(t,signal) > %stem(t(H),signal(H),'r') > %hold off > %xlabel('time','fontsize',16) > %ylabel('Amplitude','fontsize',16) > %set(gca,'fontsize',16); > %itle('Decimated signal'); > %set(gcf,'color','w'); > %set(gca,'linewidth',3) > > > %Writing code for DFT ( Discrete Fourier Transform) > %There are no weights > > N1=length(t1); > X1 = zeros(size(y1)); > > M=N1/2; > for k=1:M > freq = (-1j*2*pi*(k-1))/N1; > Wnk = exp(freq.*(0:length(t1)-1)); > X1(k)=sum(y1.*Wnk); > end > > if mod(M,2)==0 > k_dft=-N1/2:N1/2-1; % N even > else > k_dft=-(M-1)/2:(M-1)/2; % N odd > end > > > > frqaxs_dft=1/dt/N*k_dft; > > figure(4) > > plot(frqaxs_dft,fftshift(abs(X1))) > xlabel('frequency','fontsize',16) > ylabel('Power','fontsize',16) > set(gca,'fontsize',16); > title('Energy of non uniform sin signal (Positive)'); > set(gcf,'color','w'); > > > %plotting the spectrum for DFT with Weights ! > > N2=length(y1) > M2=N2/2; > if mod(M2,2)==0 > k_NDFT=-M2/2:M2/2-1; > else > k_NDFT=-(M2-1)/2:(M2-1)/2 > end > > freq_NDFT=1/dt/M2*k_NDFT; > > > > % Making the weight vector ! > > Period=max(t1)-min(t1); > deltax_n=zeros(N2,1); > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > deltax_n(N2,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > for k=2:N2-1; > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > end > > %NDFT with weights > > > N2=length(t1); > X2 = zeros(size(y1)); > M2=N2/2; > for k=1:M2 > freq = (-1j*2*pi*(k-1))/N2; > Wnk = exp(freq.*(0:length(t1)-1)); > X2(k)=sum(y1.*Wnk*deltax_n(k)); > end > > > %Contructiing a Fourier kernal > %> %Now constructing the inverse kernal > > %It should look like this finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk) * X2; > > %Although expression is like finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk+KI) * X2; > > %but author stress that if i take more number of samples than fourier coffecients than %dont need eye more it is just to suppress noise in any signal i:e for stabilization of %diagnonal > > > > > > > > result=ifft(final spectrum); > > > > > > > > > > > > > > > "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu10pv$5dp$1(a)fred.mathworks.com>... > > I read one of the paper ! It was about reconstrution of a simple sin signal ! Greag please do put your comments on it ! > > > > Objective : To reconstruct the simple sin wave through Dujhdham and Schonvillie algorithm > > > > Methodology : 1. Take DFT > > 2. Apply weights > > 3. Multiply with the inverse kernal > > 4. Simple ifft will reconstruct the samples back > > > > > > My code ( Some thing is wrong) > > > > clear all; > > close all; > > N=200; > > dt=.002; > > t=0:dt:dt*(N-1); > > signal=sin(2*pi*10*t); > > > > > > %________--------------------------------------------------- > > %_choosing percentage of decimation__ > > % Irregular interpolation > > % Randomizing the offset locations > > kn1=randperm(N); > > % Randomly eliminating some traces > > percent=0.25 > > kn2=kn1(1:ceil(percent*N)); > > % sorting random chosen traces in ascending order > > z=[sort(kn2)]; > > > > > > ------------------------------------------------------------------------- > > y1=signal; t1=t; > > > > % Finding the indices of known components > > > > y1(z)=[]; > > t1(z)=[] > > % Making the weight vector ! > > > > Period=max(t1)-min(t1); > > deltax_n=zeros(N2,1); > > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > > deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > > > for k=2:N2-1; > > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > > end > > > > %NDFT with weights > > > > > > N2=length(t1); > > X2 = zeros(size(y1)); > > > > for k=1:N2 > > freq = (-1j*2*pi*(k-1))/N2; > > Wnk = exp(freq.*(0:length(t1)-1)); > > X2(k)=sum(y1.*Wnk*deltax_n(k)); > > end > > > > % Making the weight vector ! > > > > Period=max(t1)-min(t1); > > deltax_n=zeros(N2,1); > > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > > deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > > > for k=2:N2-1; > > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > > end > > > > %DFT with weights > > > > > > N2=length(t1); > > X2 = zeros(size(y1)); > > M=N2/2; %(Author told me to always take more N and less M to make the reconstruction possible) > > for k=1:M > > freq = (-1j*2*pi*(k-1))/N2; > > Wnk = exp(freq.*(0:length(t1)-1)); > > X2(k)=sum(y1.*Wnk*deltax_n(k)); > > end > > > > %Now constructing the inverse kernal > > %It should look like this finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk) * X2; > > %Although expression is like finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk+KI) * X2; > > %but author stress that if i take more number of samples than fourier cofecients than i %dont need eye more over eye is just to suppress noise in any signal > > > > > > > > result=ifft(final spectrum); > > > > > > > > %I am doing some thing wrong here ! Author of the paper told me it should work but its not working for me ! Please guide me ! > > > > > > %If my question is not complete please let me know
From: kk KKsingh on 31 May 2010 16:09 Again a correction it should be original position instead of (0:length(t)-1)....I am sorry for this mistake ..Will keep updating Sorry for posting 4 messages ina single day "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu121b$n6o$1(a)fred.mathworks.com>... > oops another mistake inverse kernal which we multiplying with DFT spcetrum is > inv(Forward operater * W*Backward operater).... > or > inv(Forward operater * W*Backward operater+KI).... > Thanks > > "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu11m3$15t$1(a)fred.mathworks.com>... > > > I read one of the paper ! It was about reconstrution of a simple sin signal ! Greag please do put your comments on it ! > > > > > > Objective : To reconstruct the simple sin wave through Dujhdham and Schonvillie algorithm > > > > > > Methodology : 1. Take DFT > > > 2. Apply weights > > > 3. Multiply with the inverse kernal > > > 4. Simple ifft will reconstruct the samples back > > > > > > > > > My code ( Some thing is wrong) > > > > %Aki Corrected code > > clear all; > > close all; > > N=200; > > dt=.002; > > t=0:dt:dt*(N-1); > > signal=sin(2*pi*10*t); > > > > %save unisin.dat signal -ascii > > > > > > %Plotting the sin signal > > > > figure(1) > > hold on > > stem(t,signal,'r') > > plot(t,signal) > > xlabel('time','fontsize',16) > > ylabel('Amplitude','fontsize',16) > > set(gca,'fontsize',16); > > title('Sin signal'); > > set(gcf,'color','w'); > > set(gca,'linewidth',3) > > > > > > %Making Sample Missing from the Simple Sin Signal > > > > %________--------------------------------------------------- > > %_choosing percentage of decimation__ > > % Irregular interpolation > > % Randomizing the offset locations > > kn1=randperm(N); > > % Randomly eliminating some traces > > percent=0.25 > > kn2=kn1(1:ceil(percent*N)); > > % sorting random chosen traces in ascending order > > z=[sort(kn2)]; > > > > > > > > %------------------------------------------------------------------------- > > y1=signal; t1=t; > > > > % > > %h1=[1:N]; > > [%a,b]=wcommon(z,h1); > > %H=find(b-1); > > y1(z)=[]; > > t1(z)=[]; > > > > %save nonunisin.dat y1 -ascii > > %fid=fopen('nonunisin.dat','w'); > > %printf(fid,'%d\n',y1); > > %fclose(fid); > > %save nonunitime.dat t1 -ascii > > > > %figure(2) > > %hold on > > %plot(t,signal) > > %stem(t(H),signal(H),'r') > > %hold off > > %xlabel('time','fontsize',16) > > %ylabel('Amplitude','fontsize',16) > > %set(gca,'fontsize',16); > > %itle('Decimated signal'); > > %set(gcf,'color','w'); > > %set(gca,'linewidth',3) > > > > > > %Writing code for DFT ( Discrete Fourier Transform) > > %There are no weights > > > > N1=length(t1); > > X1 = zeros(size(y1)); > > > > M=N1/2; > > for k=1:M > > freq = (-1j*2*pi*(k-1))/N1; > > Wnk = exp(freq.*(0:length(t1)-1)); > > X1(k)=sum(y1.*Wnk); > > end > > > > if mod(M,2)==0 > > k_dft=-N1/2:N1/2-1; % N even > > else > > k_dft=-(M-1)/2:(M-1)/2; % N odd > > end > > > > > > > > frqaxs_dft=1/dt/N*k_dft; > > > > figure(4) > > > > plot(frqaxs_dft,fftshift(abs(X1))) > > xlabel('frequency','fontsize',16) > > ylabel('Power','fontsize',16) > > set(gca,'fontsize',16); > > title('Energy of non uniform sin signal (Positive)'); > > set(gcf,'color','w'); > > > > > > %plotting the spectrum for DFT with Weights ! > > > > N2=length(y1) > > M2=N2/2; > > if mod(M2,2)==0 > > k_NDFT=-M2/2:M2/2-1; > > else > > k_NDFT=-(M2-1)/2:(M2-1)/2 > > end > > > > freq_NDFT=1/dt/M2*k_NDFT; > > > > > > > > % Making the weight vector ! > > > > Period=max(t1)-min(t1); > > deltax_n=zeros(N2,1); > > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > > deltax_n(N2,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > > > for k=2:N2-1; > > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > > end > > > > %NDFT with weights > > > > > > N2=length(t1); > > X2 = zeros(size(y1)); > > M2=N2/2; > > for k=1:M2 > > freq = (-1j*2*pi*(k-1))/N2; > > Wnk = exp(freq.*(0:length(t1)-1)); > > X2(k)=sum(y1.*Wnk*deltax_n(k)); > > end > > > > > > %Contructiing a Fourier kernal > > %> %Now constructing the inverse kernal > > > %It should look like this finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk) * X2; > > > %Although expression is like finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk+KI) * X2; > > > %but author stress that if i take more number of samples than fourier coffecients than %dont need eye more it is just to suppress noise in any signal i:e for stabilization of %diagnonal > > > > > > > > > > > > result=ifft(final spectrum); > > > > > > > > > > > > > > > > > > > > > > > > > > > > "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu10pv$5dp$1(a)fred.mathworks.com>... > > > I read one of the paper ! It was about reconstrution of a simple sin signal ! Greag please do put your comments on it ! > > > > > > Objective : To reconstruct the simple sin wave through Dujhdham and Schonvillie algorithm > > > > > > Methodology : 1. Take DFT > > > 2. Apply weights > > > 3. Multiply with the inverse kernal > > > 4. Simple ifft will reconstruct the samples back > > > > > > > > > My code ( Some thing is wrong) > > > > > > clear all; > > > close all; > > > N=200; > > > dt=.002; > > > t=0:dt:dt*(N-1); > > > signal=sin(2*pi*10*t); > > > > > > > > > %________--------------------------------------------------- > > > %_choosing percentage of decimation__ > > > % Irregular interpolation > > > % Randomizing the offset locations > > > kn1=randperm(N); > > > % Randomly eliminating some traces > > > percent=0.25 > > > kn2=kn1(1:ceil(percent*N)); > > > % sorting random chosen traces in ascending order > > > z=[sort(kn2)]; > > > > > > > > > ------------------------------------------------------------------------- > > > y1=signal; t1=t; > > > > > > % Finding the indices of known components > > > > > > y1(z)=[]; > > > t1(z)=[] > > > % Making the weight vector ! > > > > > > Period=max(t1)-min(t1); > > > deltax_n=zeros(N2,1); > > > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > > > deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > > > > > for k=2:N2-1; > > > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > > > end > > > > > > %NDFT with weights > > > > > > > > > N2=length(t1); > > > X2 = zeros(size(y1)); > > > > > > for k=1:N2 > > > freq = (-1j*2*pi*(k-1))/N2; > > > Wnk = exp(freq.*(0:length(t1)-1)); > > > X2(k)=sum(y1.*Wnk*deltax_n(k)); > > > end > > > > > > % Making the weight vector ! > > > > > > Period=max(t1)-min(t1); > > > deltax_n=zeros(N2,1); > > > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > > > deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > > > > > for k=2:N2-1; > > > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > > > end > > > > > > %DFT with weights > > > > > > > > > N2=length(t1); > > > X2 = zeros(size(y1)); > > > M=N2/2; %(Author told me to always take more N and less M to make the reconstruction possible) > > > for k=1:M > > > freq = (-1j*2*pi*(k-1))/N2; > > > Wnk = exp(freq.*(0:length(t1)-1)); > > > X2(k)=sum(y1.*Wnk*deltax_n(k)); > > > end > > > > > > %Now constructing the inverse kernal > > > %It should look like this finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk) * X2; > > > %Although expression is like finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk+KI) * X2; > > > %but author stress that if i take more number of samples than fourier cofecients than i %dont need eye more over eye is just to suppress noise in any signal > > > > > > > > > > > > result=ifft(final spectrum); > > > > > > > > > > > > %I am doing some thing wrong here ! Author of the paper told me it should work but its not working for me ! Please guide me ! > > > > > > > > > %If my question is not complete please let me know
From: kk KKsingh on 31 May 2010 16:24
"kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu150v$2t7$1(a)fred.mathworks.com>... > Again a correction it should be original position instead of (0:length(t)-1)....I am sorry for this mistake ..Will keep updating Sorry for posting 4 messages ina single day > > "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu121b$n6o$1(a)fred.mathworks.com>... > > oops another mistake inverse kernal which we multiplying with DFT spcetrum is > > inv(Forward operater * W*Backward operater).... > > or > > inv(Forward operater * W*Backward operater+KI).... > > Thanks > > > > "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu11m3$15t$1(a)fred.mathworks.com>... > > > > I read one of the paper ! It was about reconstrution of a simple sin signal ! Greag please do put your comments on it ! > > > > > > > > Objective : To reconstruct the simple sin wave through Dujhdham and Schonvillie algorithm > > > > > > > > Methodology : 1. Take DFT > > > > 2. Apply weights > > > > 3. Multiply with the inverse kernal > > > > 4. Simple ifft will reconstruct the samples back > > > > > > > > > > > > My code ( Some thing is wrong) > > > > > > %Aki Corrected code > > > clear all; > > > close all; > > > N=200; > > > dt=.002; > > > t=0:dt:dt*(N-1); > > > signal=sin(2*pi*10*t); > > > > > > %save unisin.dat signal -ascii > > > > > > > > > %Plotting the sin signal > > > > > > figure(1) > > > hold on > > > stem(t,signal,'r') > > > plot(t,signal) > > > xlabel('time','fontsize',16) > > > ylabel('Amplitude','fontsize',16) > > > set(gca,'fontsize',16); > > > title('Sin signal'); > > > set(gcf,'color','w'); > > > set(gca,'linewidth',3) > > > > > > > > > %Making Sample Missing from the Simple Sin Signal > > > > > > %________--------------------------------------------------- > > > %_choosing percentage of decimation__ > > > % Irregular interpolation > > > % Randomizing the offset locations > > > kn1=randperm(N); > > > % Randomly eliminating some traces > > > percent=0.25 > > > kn2=kn1(1:ceil(percent*N)); > > > % sorting random chosen traces in ascending order > > > z=[sort(kn2)]; > > > > > > > > > > > > %------------------------------------------------------------------------- > > > y1=signal; t1=t; > > > > > > % > > > %h1=[1:N]; > > > [%a,b]=wcommon(z,h1); > > > %H=find(b-1); > > > y1(z)=[]; > > > t1(z)=[]; > > > > > > %save nonunisin.dat y1 -ascii > > > %fid=fopen('nonunisin.dat','w'); > > > %printf(fid,'%d\n',y1); > > > %fclose(fid); > > > %save nonunitime.dat t1 -ascii > > > > > > %figure(2) > > > %hold on > > > %plot(t,signal) > > > %stem(t(H),signal(H),'r') > > > %hold off > > > %xlabel('time','fontsize',16) > > > %ylabel('Amplitude','fontsize',16) > > > %set(gca,'fontsize',16); > > > %itle('Decimated signal'); > > > %set(gcf,'color','w'); > > > %set(gca,'linewidth',3) > > > > > > > > > %Writing code for DFT ( Discrete Fourier Transform) > > > %There are no weights > > > > > > N1=length(t1); > > > X1 = zeros(size(y1)); > > > > > > M=N1/2; > > > for k=1:M > > > freq = (-1j*2*pi*(k-1))/N1; > > > Wnk = exp(freq.*(0:length(t1)-1)); > > > X1(k)=sum(y1.*Wnk); > > > end > > > > > > if mod(M,2)==0 > > > k_dft=-N1/2:N1/2-1; % N even > > > else > > > k_dft=-(M-1)/2:(M-1)/2; % N odd > > > end > > > > > > > > > > > > frqaxs_dft=1/dt/N*k_dft; > > > > > > figure(4) > > > > > > plot(frqaxs_dft,fftshift(abs(X1))) > > > xlabel('frequency','fontsize',16) > > > ylabel('Power','fontsize',16) > > > set(gca,'fontsize',16); > > > title('Energy of non uniform sin signal (Positive)'); > > > set(gcf,'color','w'); > > > > > > > > > %plotting the spectrum for DFT with Weights ! > > > > > > N2=length(y1) > > > M2=N2/2; > > > if mod(M2,2)==0 > > > k_NDFT=-M2/2:M2/2-1; > > > else > > > k_NDFT=-(M2-1)/2:(M2-1)/2 > > > end > > > > > > freq_NDFT=1/dt/M2*k_NDFT; > > > > > > > > > > > > % Making the weight vector ! > > > > > > Period=max(t1)-min(t1); > > > deltax_n=zeros(N2,1); > > > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > > > deltax_n(N2,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > > > > > for k=2:N2-1; > > > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > > > end > > > > > > %NDFT with weights > > > > > > > > > N2=length(t1); > > > X2 = zeros(size(y1)); > > > M2=N2/2; > > > for k=1:M2 > > > freq = (-1j*2*pi*(k-1))/N2; > > > Wnk = exp(freq.*(0:length(t1)-1)); > > > X2(k)=sum(y1.*Wnk*deltax_n(k)); > > > end > > > > > > > > > %Contructiing a Fourier kernal > > > %> %Now constructing the inverse kernal > > > > %It should look like this finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk) * X2; > > > > %Although expression is like finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk+KI) * X2; > > > > %but author stress that if i take more number of samples than fourier coffecients than %dont need eye more it is just to suppress noise in any signal i:e for stabilization of %diagnonal > > > > > > > > > > > > > > > > result=ifft(final spectrum); > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > "kk KKsingh" <akikumar1983(a)gmail.com> wrote in message <hu10pv$5dp$1(a)fred.mathworks.com>... > > > > I read one of the paper ! It was about reconstrution of a simple sin signal ! Greag please do put your comments on it ! > > > > > > > > Objective : To reconstruct the simple sin wave through Dujhdham and Schonvillie algorithm > > > > > > > > Methodology : 1. Take DFT > > > > 2. Apply weights > > > > 3. Multiply with the inverse kernal > > > > 4. Simple ifft will reconstruct the samples back > > > > > > > > > > > > My code ( Some thing is wrong) > > > > > > > > clear all; > > > > close all; > > > > N=200; > > > > dt=.002; > > > > t=0:dt:dt*(N-1); > > > > signal=sin(2*pi*10*t); > > > > > > > > > > > > %________--------------------------------------------------- > > > > %_choosing percentage of decimation__ > > > > % Irregular interpolation > > > > % Randomizing the offset locations > > > > kn1=randperm(N); > > > > % Randomly eliminating some traces > > > > percent=0.25 > > > > kn2=kn1(1:ceil(percent*N)); > > > > % sorting random chosen traces in ascending order > > > > z=[sort(kn2)]; > > > > > > > > > > > > ------------------------------------------------------------------------- > > > > y1=signal; t1=t; > > > > > > > > % Finding the indices of known components > > > > > > > > y1(z)=[]; > > > > t1(z)=[] > > > > % Making the weight vector ! > > > > > > > > Period=max(t1)-min(t1); > > > > deltax_n=zeros(N2,1); > > > > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > > > > deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > > > > > > > for k=2:N2-1; > > > > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > > > > end > > > > > > > > %NDFT with weights > > > > > > > > > > > > N2=length(t1); > > > > X2 = zeros(size(y1)); > > > > > > > > for k=1:N2 > > > > freq = (-1j*2*pi*(k-1))/N2; > > > > Wnk = exp(freq.*(0:length(t1)-1)); > > > > X2(k)=sum(y1.*Wnk*deltax_n(k)); > > > > end > > > > > > > > % Making the weight vector ! > > > > > > > > Period=max(t1)-min(t1); > > > > deltax_n=zeros(N2,1); > > > > deltax_n(1,:)=(t1(2)-(t1(N2)-Period))/2; > > > > deltax_n(N1,:)=((t1(1)+Period)-(t1(N2-1)))/2; > > > > > > > > for k=2:N2-1; > > > > deltax_n(k,:)= ((t1(k+1)-t1(k-1))/2); > > > > end > > > > > > > > %DFT with weights > > > > > > > > > > > > N2=length(t1); > > > > X2 = zeros(size(y1)); > > > > M=N2/2; %(Author told me to always take more N and less M to make the reconstruction possible) > > > > for k=1:M > > > > freq = (-1j*2*pi*(k-1))/N2; > > > > Wnk = exp(freq.*(0:length(t1)-1)); > > > > X2(k)=sum(y1.*Wnk*deltax_n(k)); > > > > end > > > > > > > > %Now constructing the inverse kernal > > > > %It should look like this finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk) * X2; > > > > %Although expression is like finalspectrum= inv(Wnk'*diag(deltax_n)*Wnk+KI) * X2; > > > > %but author stress that if i take more number of samples than fourier cofecients than i %dont need eye more over eye is just to suppress noise in any signal > > > > > > > > > > > > > > > > result=ifft(final spectrum); > > > > > > > > > > > > > > > > %I am doing some thing wrong here ! Author of the paper told me it should work but its not working for me ! Please guide me ! > > > > > > > > > > > > %If my question is not complete please let me know Hello ! I think i need to leave my habit of top posting :(! Any ways can some have a look on my code and suggest some correction ! My ultimate aim right now is to esimate less number of fourier cofficients than number of samples and i am doing totally wrong some where ! So it will be great if some body can suggest me correction |