matlab code to c

13 views (last 30 days)
Nandini
Nandini on 22 Jun 2022
Answered: Ishaan Mehta on 27 Jun 2022
please can anyone help me to solve this ?? m hving a matlab code convert this code into C programm
Decoding VMD algorithm
Step01: Initialization of all required parameters based on the length of the signal
Step02: Reduce edge effect by mirroring signal;
Step03: Signal frequency domain with full bandwidth
Step04: Get half of the bandwidth
Step04: FFT for initial IMFs and get half of bandwidthl;
Step05: Obtain Frequency vector
Step06: Get the initial central frequencies
Step07: Optimization iterations
Step08: Convert to time domain signal
Step09: Calculate residual
%%--------------------- Step01: --------------------------------
PenaltyFactor=100;
LMUpdateRate=0.0100;
AbsoluteTolerance=5.0000e-06;
RelativeTolerance=0.0050;
MaxIterations=1000;
InitialIMFs=zeros(length(x),5);
InitialLM=zeros(length(x)+1,1);
CentralFrequencies= [];
InitializeMethod='peaks';
FFTLength=2*length(x);
NumIMFs= 5;
SignalLength= length(x);
HalfSignalLength=length(x)/2;
MirroredSignalLength=2*length(x);
DataType= 'double';
NumHalfFreqSamples= length(x)+1;
Display= 0;
%%
nfft = FFTLength;
penaltyFactor = PenaltyFactor;
numIMFs = NumIMFs;
relativeDiff = inf;
absoluteDiff = relativeDiff;
tau = LMUpdateRate; % Lagrange multiplier update rate
%%--------------------- Step02 & Step 03-------------------------
xr = [x(HalfSignalLength:-1:1); x; x(SignalLength:-1:ceil(SignalLength/2)+1)];
y = fft(xr,FFTLength);
sigFDFull = y;
% Get half of the bandwidth
sigFD = sigFDFull(1:NumHalfFreqSamples);
%%--------------------- Step 04 --------------------------------
initIMFfdFull = fft(InitialIMFs,nfft);
initIMFfd = initIMFfdFull(1:NumHalfFreqSamples,:) + eps;
IMFfd = initIMFfd;
sumIMF = sum(IMFfd,2);
LM = InitialLM(:); % Lagrange Multiplier
%% Frequency vector from [0,0.5) for odd nfft and [0,0.5] for even nfft
%%--------------------- Step 05 --------------------------------
f = (0:(nfft/2))/nfft;
%%--------------------- Step 06 --------------------------------
%% Get the initial central frequencies
x=abs(sigFD);
BW = 2/FFTLength; % bandwidth of signal
minBWGapIndex = 2*BW/f(2);
x(x<mean(x)) = mean(x);
TF = islocalmax(x,'MinSeparation',minBWGapIndex);
pkst = x(TF);
locst = f(TF);
numpPeaks = length(pkst);
% Check for DC component
if x(1) >= x(2)
pks = zeros(numpPeaks+1,1);
locs = pks;
pks(2:length(pkst)+1) = pkst;
locs(2:length(pkst)+1) = locst;
pks(1) = x(1);
locs(1) = f(1);
else
pks = zeros(numpPeaks,1);
locs = pks;
pks(1:length(pkst)) = pkst;
locs(1:length(pkst)) = locst;
end
[~,index] = sort(pks,'descend');
centralFreq = 0.5*rand(NumIMFs,1);
% Check if the number of peaks is less than number of IMFs
if length(locs) < NumIMFs
centralFreq(1:length(locs(index))) = locs;
else
centralFreq(1:NumIMFs) = locs(index(1:NumIMFs));
end
%%
iter = 0;
f=f';
initIMFNorm = abs(initIMFfd).^2;
normIMF = zeros(size(initIMFfd,1),size(initIMFfd,2));
%%--------------------- Step 07 --------------------------------
while (iter < MaxIterations && (relativeDiff > RelativeTolerance ||...
absoluteDiff > AbsoluteTolerance))
for kk = 1:numIMFs
sumIMF = sumIMF - IMFfd(:,kk);
IMFfd(:,kk) = (sigFD - sumIMF + LM/2)./...
(1+penaltyFactor*(f - centralFreq(kk)).^2);
normIMF(:,kk) = abs(IMFfd(:,kk)).^2;
centralFreq(kk) = (f.'*normIMF(:,kk))/sum(normIMF(:,kk));
sumIMF = sumIMF + IMFfd(:,kk);
end
LM = LM + tau*(sigFD-sumIMF);
absDiff = mean(abs(IMFfd-initIMFfd).^2);
absoluteDiff = sum(absDiff);
relativeDiff = sum(absDiff./mean(initIMFNorm));
% Sort IMF and central frequecies in descend order
% In ADMM, the IMF with greater power will be substracted first
[~,sortedIndex] = sort(sum(abs(IMFfd).^2),'descend');
IMFfd = IMFfd(:,sortedIndex);
centralFreq = centralFreq(sortedIndex(1:length(centralFreq)));
initIMFfd = IMFfd;
initIMFNorm = normIMF;
iter = iter + 1;
end
%%--------------------- Step 08 --------------------------------
%% Convert to time domain signal
% Transform to time domain
IMFfdFull = complex(zeros(nfft,numIMFs));
IMFfdFull(1:size(IMFfd,1),:) = IMFfd;
if ~mod(FFTLength,2)
IMFfdFull(size(IMFfd,1)+1:end,:) = conj(IMFfd(end-1:-1:2,:));
else
IMFfdFull(size(IMFfd,1)+1:end,:) = conj(IMFfd(end:-1:2,:));
end
[~,index] = sort(centralFreq,'descend');
%%
z=IMFfdFull(:,index);
xr = real(ifft(z,FFTLength));
IMFs_without_inbuild = xr(HalfSignalLength+1:MirroredSignalLength-HalfSignalLength,:);
%%--------------------- Step 09 --------------------------------
residual_without_inbuild = PPGblr1 - sum(IMFs_without_inbuild,2);

Answers (1)

Ishaan Mehta
Ishaan Mehta on 27 Jun 2022
Hi Nandini
I understand that you want to convert your existing MATLAB code to C code.
I suggest you to refer to Generate C Code by Using the MATLAB Coder App article for the same.
Hope this helps
Ishaan Mehta

Categories

Find more on Get Started with Signal Processing Toolbox in Help Center and File Exchange

Tags

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!