Need help to find true Signal to noise ratio (SNR) for acceleration signal

Question:

I am quite new to matlab, and I have gotten the task to find the true Signal to noise ratio for the original acceleration signal. For the signal to noise ratio (SNR1) for the noisy acceleration signal, we got a value of 6.93.

So what I am struggling with is how to find the sum harmonics I need to use the code I was thinking. Does anybody know how to do that? or is it another better way?

file_path=uigetdir;
cd(file_path);
files=dir("acc_data*.mat");
load Index_table.mat

Fig=1
Fs=100; %sampling frequency
F= Fs/2;

for i=1:length(files)
    load(files(i).name);%task 2
    indx=find(time>=Index_tabel2{i,1}& time<=Index_tabel2{i,2}); %task 3
    N=size(acc,2);
    Win=10; %why 10 as a window size?
    win_center=round(Win/2);
    PSD_acc=[];Freq_acc=[];PSD_acc1=[];Freq_acc1=[];
    
    for n=1:N
        sd_acc=std(acc(indx,n)); %Standard deviation of acc indx
        m_acc=mean(acc(indx,:)); %mean ...
        acc_noise= (0.3*sd_acc)*randn(size(indx,1),n); %task 4:creating noisy version of signal
        
        NOISYsignal=acc_noise+acc(indx,n); %Add noise with amplitude 0.3*std(acc)to each original acceleration signal 
        NOISYsignal_sub=NOISYsignal-mean(NOISYsignal);
        ORIGINALacc_sub=acc(indx,:)-m_acc;
        
        [PSD,Freq]=pwelch(ORIGINALacc_sub,1000,[],F,Fs); %PSD of original signal
        [PSD1,Freq1]= pwelch(NOISYsignal_sub,1000,[],F,Fs);%PSD of noisy+original signal
        
        PSD_acc=[PSD_acc,PSD]; 
        PSD_acc1=[PSD_acc1,PSD1];
        Freq_acc=[Freq_acc,Freq]; 
        Freq_acc1=[Freq_acc1,Freq1];
        
        
        Cutoff = 20;
        Indx_signal = find(Freq1<=Cutoff);
        Indx_noise = find(Freq1>Cutoff);
        PSD_signal1 = sum(PSD1(Indx_signal));
        PSD_noise1 = sum(PSD1(Indx_noise));
        SNR1 = PSD_signal1/PSD_noise1;
        
        %True SNR (this is the code we will use to find true SNR
        %     true_noise = signal-sum_harmonics;
        %     var_signal_true = var(sum_harmonics);
        %     var_noise_true = var(signal-sum_harmonics);
        %     SNR_true = var_signal_true/var_noise_true;
        
        
        %moving average
        for m=1:length(NOISYsignal_sub)-Win
            smooth_acc(m)=mean(NOISYsignal_sub(m:m+Win));
            win_time(m)=time(indx(m+win_center));
            
            
        end  
        
        
    end
    if Fig==1
        
        figure;%task 8: moving average
        plot(time(indx),NOISYsignal_sub,'b-');hold on %plotting of noisy signal: by taking minus mean--> zero
        plot(time(indx),acc(indx),'r-',"LineWidth",1); %plotting of original data: there is an off set in the original data, it happens due to incorrect offset
        plot(win_time,smooth_acc,'g',"LineWidth",2); %plotting of moving average in green
        xlabel('time');ylabel('acceleration');
        title('Noisy acceleration with moving average')
        
        figure;%task 5b
        for n=1:N
            subplot(2,2,n)
            plot(Freq_acc(:,n),PSD_acc(:,n));hold on
            plot(Freq_acc1(:,n),PSD_acc1(:,n));
            title(Labels_acc{n});
            xlabel('Frequency(Hz)');ylabel('PSD');
            legend('Original Signal','Noisy signal');
        end
        
    end
    
end

[hanna Eid]

How do you make PSD frequency into a signal?

Answer 1

[Abhimenyu]

Hi Hanna,

I understand that you want to find the true Signal-to-noise ratio of the acceleration signal using the sum of harmonics. The “harmonic” function of MATLAB can be used to find the sum of harmonics. The “harmonic” function is defined as the sum of reciprocals of the first “n” positive integers where “n” is the input.

To find the sum of harmonics of a signal, first MATLAB’s “fft” function can be used to compute the discrete Fourier transform of the signal, and then the “harmonic” function can be used to sum the reciprocals of the magnitudes of the frequency components. Please refer to the example code below to understand more:

% Compute the discrete Fourier transform of the signal 
X = fft(signal); 

% Get the magnitudes of the frequency components 
Mag = abs(X); 

% Sum the reciprocals of the magnitudes using the "harmonic" function 
sum_harmonics = harmonic(M); 

The variable “sum_harmonics” will contain the sum of harmonics of the signal. This value can be used to compute the true signal-to-noise ratio (SNR) of the signal according to the code mentioned by you.

Please refer to the below given MATLAB documentation links to explore more on the “harmonic” and the “fft” function respectively:

• https://www.mathworks.com/help/symbolic/sym.harmonic.html
• https://www.mathworks.com/help/matlab/ref/fft.html

I hope this helps to resolve the query.

Thanks,

Abhimenyu