How do you make PSD frequency into a signal?

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

8 views (last 30 days)

Show older comments

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

### Answers (1)

Abhimenyu
on 8 Dec 2023

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

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!