Creating a Gaussian Filter Transformation Matrix for 1D Signals

Question

I have a signal of size 120x1. I am trying to create a 1D Gaussian filter (G) such that when it is multiplied by the signal it generates the filtered signal.

sFilt = G*s

I know that MATLAB filter function can apply this transformation by:

i = -30:30;
sigma = 2;
gaussFilter = exp(- i.^ 2 / (2 * sigma ^ 2));
gaussFilter = gaussFilter / sum (gaussFilter); % normalize
s = rand(120,1); % Example signal
sFilt = filter(gaussFilter,1, s); % Example filtered signal

However, I want to find the kernel matrix by myself rather than MATLAB internally finds and applies it to my signal.

Here is my signal for which I want to find G.

s = [...
       0.00694444000000000  0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.930556000000000   3.38767207896461    2.96157398017693    -0.333333000000000  1.60762896805897    2.00694000000000    2.00694000000000    2.00694000000000    1.28472000000000    -1.43744631449631   0.944444000000000   2.59392295081968    2.99306000000000    2.99306000000000    2.99306000000000    2.99306000000000    2.99306000000000    2.99306000000000    2.99306000000000    2.02922488524590    -0.812570065573777  -1.18993444717445   -0.782122839475842  0.640152000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   2.01873464373465    4.49027190008191    5.18750000000000    5.18750000000000    5.18750000000000    5.18750000000000    5.18750000000000    4.59722000000000    2.59722000000000    3.97561321867322    3.21314756756758    0.604167000000000   0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.111111000000000   0.458646085995084   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000
    ]';

Here is the filtered signal with added noise.

sFilt = [...
       0.0138889000000000   0.00208355909090909 -0.00329849545454587    0.201389000000000   0.861111500000000   1.86806000000000    2.44947936559843    2.52500275098103    2.15972000000000    2.04861000000000    2.14807844644317    2.16668303352412    1.78457304987735    1.27778000000000    1.00173350000000    1.57634828431373    2.33333000000000    3.21248677841374    3.55556000000000    3.64930500000000    3.71804717320261    3.81944000000000    3.49947244281046    2.76389000000000    1.65451500000000    0.470087078495503   -0.229167000000000  -0.208333000000000  0.354167000000000   0.789952043336059   1.08820161079313    1.19444000000000    1.19444000000000    1.07082559280458    1.09288587898610    1.19444000000000    1.20139000000000    1.23611000000000    1.45729755519215    2.02778000000000    3.14999395748160    4.65278000000000    5.78193152902698    6.30624161896975    6.45139000000000    6.34444931316435    5.92361000000000    5.44723022076860    4.78127380212592    4.18750000000000    3.06944000000000    1.72222000000000    0.791893542109568   0.312170886345054   0.0555556000000000  0.00694444000000000 0.0482593386437905  0.0524382556009807  -0.0647580647058824 -0.0694444000000000 0.0486111000000000  0.0249897492232223  0.0356016021241839  0.256944000000000   0.361111000000000   0.573964112019624   0.681776461161079   0.674996493049878   0.583855436631235   0.625000000000000   0.701389000000000   0.701389000000000   0.673611000000000   0.729167000000000   0.701389000000000   0.659722000000000   0.583333000000000   0.645833000000000   0.687500000000000   0.701389000000000   0.625000000000000   0.625000000000000   0.701389000000000   0.701389000000000   0.629343798855273   0.583333000000000   0.625000000000000   0.645833000000000   0.631944000000000   0.460070669664759   0.665622201962387   0.527778000000000   0.562500000000000   0.520833000000000   0.562500000000000   0.541667000000000   0.645833000000000   0.708333000000000   0.581075065359477   0.579855317512275   0.548611000000000   0.666667000000000   0.548611000000000   0.618056000000000   0.576389000000000   0.590278000000000   0.645833000000000   0.666667000000000   0.554531401639344   0.548611000000000   0.576389000000000   0.583333000000000   0.656560276422764   0.584357622950820   0.559938565573770   0.465278000000000   0.583333000000000   0.659722000000000   0.687500000000000   0.673611000000000
      ]';

I tried the following code to obtain G but it does not generate the correct filtered signal:

G = zeros(180); % Zero padding G
i = -30:30; 
sigma = 2;
gaussFilter = exp(- i.^ 2 / (2 * sigma ^ 2));
gaussFilter = gaussFilter / sum (gaussFilter); % normalize
for i = 31 : 150
    G(i, i-30:i+30) = gaussFilter;
end
G = G(31:150,31:150); % Truncate to size 120x120

As can be seen, G*s is different from sFilt regardless of the noise. My goal is to be able to find G such that G*s = sFilt regardless of the noise and the filtered signal has approximately the same amplitude.

May I also know if the way I created G is correct? If not, how should I build G for the signal s? Given s and sFilt, is it correct to find G = sFilt * pinv(s), which I am afraid of getting a singular matrix?

I obtain exactly the same result. Are you sure you used the same input s for both? Maybe a single minimal reproducible example may clarify your problem (including the plotting code). — m7913d
s is used in G*s to compare with sFilt. I don't get the correct amplitude when I apply G*s. — user578
I can't reproduce your graph. I'm getting exactly the same between them both. BTW, I had to change s = rand(500,1); to s = rand(120,1); to get the multiplication to work. — rayryeng
I just want to know how to create G. Those codes are just examples to clarify my question. — user578

m7913d m7913d · Accepted Answer · 2017-09-07T06:40:21

I think you should have made a typo, as I obtain exactly the same result for both methods:

s = [...
       0.00694444000000000  0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.930556000000000   3.38767207896461    2.96157398017693    -0.333333000000000  1.60762896805897    2.00694000000000    2.00694000000000    2.00694000000000    1.28472000000000    -1.43744631449631   0.944444000000000   2.59392295081968    2.99306000000000    2.99306000000000    2.99306000000000    2.99306000000000    2.99306000000000    2.99306000000000    2.99306000000000    2.02922488524590    -0.812570065573777  -1.18993444717445   -0.782122839475842  0.640152000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   0.909722000000000   2.01873464373465    4.49027190008191    5.18750000000000    5.18750000000000    5.18750000000000    5.18750000000000    5.18750000000000    4.59722000000000    2.59722000000000    3.97561321867322    3.21314756756758    0.604167000000000   0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.00694444000000000 0.111111000000000   0.458646085995084   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000   0.506944000000000
    ];

i = -30:30;
sigma = 2;
gaussFilter = exp(- i.^ 2 / (2 * sigma ^ 2));
gaussFilter = gaussFilter / sum (gaussFilter); % normalize
sFilt = filter(gaussFilter,1, s); % Example filtered signal

G = zeros(180); % Zero padding G
for i = 31 : 150
    G(i, i-30:i+30) = gaussFilter;
end
G = G(31:150,31:150); % Truncate to size 120x120

figure
hold on
plot(s)
plot(sFilt(31:end))
plot(G*s')

Creating a Gaussian Filter Transformation Matrix for 1D Signals

1 Answers