Inverse MNF transform

[AdamVarley30]

Edzer,
I know that this should be basic knowledge, so I apologise in advance.

I was contacting you to inquire how you can rebuild/rescale MNF results to reconstruct a clearer image/ denoise data. Presumably this is very similar to a conventional PCA transform (see below). However, I cannot for the life of me get it to work.

t(t(pca$x[,1:20] %% t(pca$rotation[,1:20])) pca$scale + pca$center)

Any help would be fantastic

Kind regards

Adam Varley

[edzer]

If this is still relevant, please provide a reproducible example.

[AdamVarley30]

Hello Edzer,

Still very relevant!!!

Thanks for getting back. I have a number of noisy gamma-ray spectra (a few thousand or so) and i want to reconstruct or denoise each measurement using the entire dataset employing MNF. This has been performed by a number of other studies in the past using 'unknown' software (probably ENVI).

Essentially, I have a matrix with energy channels in the columns and individual measurements in each row. I want use your MNF algorithm to best estimate the signal component (lets assume it lies within the top 40 components). I then want to reproduce each spectrum using only these 40 components using the corresponding eigen values for each observation, thus discarding of some of the some of the counting noise that is distributed across each spectrum within the lower order components. I hope that make sense?

As I say, i can do it using conventional PCA (using the code i provided a year ago), but I'm unsure of how to do it using your algorithm, It should be theoretically possible, I just don't have the mathematical prowess to be able to do it.

I can provide some code if you need an example?

Thanks again for getting back

[AdamVarley30]

Hi Edzer,

It would be great if you could help me.

Thanks

Adam

[edzer]

Yes, I asked for a reproducible example, something simple, with data and code.

[AdamVarley30]

Hey Edzer here is my example. I am really lost on reconstructing a clean dataset.

Below is the code and i have attached a .txt with the data. So should be set up.
Any help would be enormous!!

Thanks again for your time I know you are a busy man.

p.s. on a side note the sf package has revolutionized R!!! Well done sir, you deserve a knighthood.

forEdzer.txt

rm(list=ls())
setwd("C:/")
Field.data <- read.table("forEdzer.txt",sep = " ")

Isolate useful energy channels

startwindow <- 50
Field.data.1 <- Field.data[,startwindow:1010]

Base PCA

pca <- prcomp(Field.data.1,scale. = T)

use the first 12 PC to reconstruct each spectra

transposed <- t(t(pca$x[,1:12] %% t(pca$rotation[,1:12])) pca$scale + pca$center)

Plot of noisy vs "cleaned"

spec.no <- 10
noisy <- as.numeric(Field.data.1[spec.no,])
cleaned <- as.numeric(transposed[spec.no,])
plot(noisy,log="y",ylim=c(0.05,20),las=1,font.axis=2,ylab="counts per second",xlab = "Channel number")
lines(cleaned)

Using spacetime

library(spacetime)

Base PCA

mnf_result <- mnf(x = as.matrix(Field.data.1,))

use the first 12 PC from MNF to reconstruct each spectra

transposed_mnf <- t(t(mnf_result$x[,1:12] %*% t(mnf_result$rotation[,1:12])))

spec.no <- 10
noisy <- as.numeric(Field.data.1[spec.no,])
cleaned <- as.numeric(transposed_mnf[spec.no,])

The results look terrible (i think it is something to do with the scaling!!)

plot(noisy,log="y",ylim=c(0.05,20),las=1,font.axis=2,ylab="counts per second",xlab = "Channel number")
lines(cleaned)

[AdamVarley30]

Edzer, is the above code reproducible?