4
$\begingroup$

I'm trying to simulate the 2nu2 band of HCN in MATLAB. I generate a stick spectrum from the data I have (taken from the HITRAN database) and then use conv() to apply my line shape function. However, when I plot the data, the broadened spectrum does not line up with the stick spectrum at all. I cannot figure out what it is that I'm doing wrong or what I need to change to align the spectra. Does anyone know what is I need to fix or a resource I could use to solve this issue?

This is the spectrum I've simulated:

Simulated HCN spectra

and the function I use to generate the spectrum is

function f = func()
HCN2nu2 = importdata('I:/HCN 2nu2.txt');
HCN = HCN2nu2(1:111,1:2);
AX = HCN2nu2(1:111,2);
v = [];
s = [];
L = [];
q = 1;
p = HCN(q);
increase = false;
for i = 6200:0.01:6700
    for x = 1:111
        if i == round(HCN(x,1), 2)
            z = AX(x);
            increase = true;
            break
        else
            z = 0;
        end
    end
    v = [v; i];
    s = [s; z];
    if increase == true
        q = q + 1;
        increase = false;
        if q > 111
            q = 111;
        end
    end
    t = (1/(1 + ((i - p)/0.5)^2));
    L = [L; t];
end
hold on
plot(v,s)
k = conv(L,s, 'same');
plot(v,k)
ylim([0 8e-21]);
hold off
end

Any help would be greatly appreciated

Link to the data file used to generate this spectrum: https://github.com/tkh97/HCN-2nu2.git

$\endgroup$
  • 2
    $\begingroup$ This is arguably off-topic for Chemistry.SE. Although the code is chemistry-related, it isn't actually asking about chemistry, it's a programming exercise. My suggestion would be to ask this on MATLAB Answers. There are lots of qualified people there. And you should probably make that txt file available, or else it's rather difficult to figure out what's going on. I know enough MATLAB to do work with it, but I am having a difficult time following your code since I have no clue what vectors HCN and AX are. $\endgroup$ – orthocresol Jul 27 at 4:08
  • 1
    $\begingroup$ Signal Processing stack exchange is a good resource for these questions but his/her error is rather simple to fix "conceptually". It is not a Matlab issue, which is doing the convolution correctly, but rather it is a minor erroneous math. $\endgroup$ – M. Farooq Jul 27 at 4:29
4
$\begingroup$

Update : In the given data by the OP, the x scale is not evenly spaced so the discrete convolution approach will not work to generate a spectrum until and unless your data is evenly spaced! The trick mentioned by EdV is a very nice shortcut. Here, if you really want to play around with convolution via MATLAB, then read below.

Solution: You can satisfy yourself as follows:

  1. Make sure you have a even spacing of the x-axis. Your uneven axis threw me off, and it took a while to discover why it is not working. I wish you mentioned it earlier.
  2. Your sticks are basically scaled delta functions as Prof. Ed said.
  3. Generate a vector for a single Lorentzian peak, whose center is zero, which is also evenly spaced. The x-axis of the Lorentzian would be symmetric [-length (wavenumber/2):Sampling rate:length(wavenumber/2)-Sampling rate].
  4. Perform convolution of your spectrum and this Lorentzian.

Let us start with a Dirac delta centered at point $x.$ I don't know how to draw a displace Dirac delta in MATLAB. I manually made it in an Excel file. If you convolute it with a zero centered Lorentzian, then you will get a Lorentzian centered at $x.$ This is called the sifting property.

Basically, all you have to ensure is that your Lorentzian is centered at zero. First, play with convolution using single stick. Once you correct the code, so that the line position does not change, apply it on the entire spectrum.

Convolution

A great book is Bracewell's Fourier Transform And Its Applications. It is all there. From MIT 2.14 / 2.140 Analysis and Design of Feedback Control Systems, Spring 2007 handout:

1.1 The “Sifting” Property of the Impulse

When an impulse appears in a product within an integrand, it has the property of ”sifting” out the value of the integrand at the point of its occurrence:

$$\int_{-\infty}^\infty f(t)δ(t - a)\,\mathrm dt = f(a)\tag{5}$$

This is easily seen by noting that $δ(t − a)$ is zero except at $t = a,$ and for its infinitesimal duration $f(t)$ may be considered a constant and taken outside the integral, so that

$$\int_{-\infty}^\infty f(t)δ(t - a)\,\mathrm dt = f(a)\int_{-\infty}^\infty δ(t - a)\,\mathrm dt = f(a)\tag{6}$$

from the unit area property.

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ I do suspect you are quite right; however, it's really hard for me to verify this (and thus upvote) since those vectors are coming from a black box .txt file. Very nice suggestion on Sig Proc SE though. $\endgroup$ – orthocresol Jul 27 at 4:50
  • 1
    $\begingroup$ @orthocresol, I am recently writing a paper in on a new approach towards Fourier self-deconvolution for resolving spectral lines. These shifting issues during convolution came up again and again and it took me at least a year to sort them out. OP is seeing exactly the same issue. Many thanks to Prof. Ed as well. I cannot parse his Matlab code. I am rather an elementary user of Matlab. $\endgroup$ – M. Farooq Jul 27 at 4:55
  • 1
    $\begingroup$ (+1) Funny thing about this one: the stick figure spectrum is just a scaled set of “delta functions”, and convolution with a “delta function” is the identity operation, so it looks like all that is necessary is to place a “stick height”-scaled Lorentzian (with 1 wavenumber FWHM) at each of the sticks in the raw spectrum. $\endgroup$ – Ed V Jul 27 at 11:49
2
$\begingroup$

This is an elaborated explanation of what commented upon. As stated by others, the question likely should have been posted to another stack exchange, e.g., signal processing. Furthermore, the OP really just wants to know why their Matlab code results in a shifted HCN spectrum, but is otherwise correct. The answer by M. Farooq, which I upvoted, provides the OP with what they need to consider as they revise and test their program. Here, however, I simply show how to do solve the problem without doing explicit convolution.

To begin, consider a unit height Lorentzian line profile function:

$$Lor(x;L,W) = \frac{1}{1+4[(x-L)/W]^2} \tag{1}$$

where x is the independent variable, L is the peak's centroid location and W is the FWHM (full width at half maximum height) of the peak. Units: x, L and W are in wavenumbers, i.e., $cm^{-1}$. From the legend of the OP's figure, W = 1 $cm^{-1}$, so equation (1) reduces to

$$Lor(x;L,1) = \frac{1}{1+4(x-L)^2} \tag{2}$$

If L = 0, then the Lorentzian peak is zero-centered:

$$Lor(x;0,1) = \frac{1}{1+4x^2} \tag{3}$$

The data file linked by the OP gives 111 L,h pairs. These are the respective locations and heights of the 111 "sticks" in the OP's stick figure HCN spectrum. Plotting these yields the following result:

HCN stick spectrum

Now for the short cut. Each of the "sticks" in the HCN spectrum may be considered as a scaled and shifted "delta function". Thus the $i^{th}$ "stick", having $L_i$ location and $h_i$ height, is

$$stick_i = h_i \times \delta(x-L_i) \tag{4}$$

But convolution with a "delta function" is the identity operation, so convolution of a zero-centered Lorentzian peak profile, with a scaled and shifted "delta function", simply results in an $h_i$-scaled Lorentzian peak at location $L_i$. Hence

$$\frac{h_i}{1+4(x-L_i)^2} = Lor(x;0,1) * h_i \times \delta(x-L_i) \tag{5}$$

where convolution is denoted by "*".

Now all that is necessary is to evaluate 111 Lorentzians, one per $L_i,h_i$ pair, as per the left-hand side of equation (5), and then add them all up. I did this with a simple computer program (not Matlab), with x ranging from 6200.00 to 6700.00, in increments of 0.01 $cm^{-1}$. The resulting spectrum is

HCN Lorentzians spectrum

Even though the $h_i$ values are unevenly spaced and given to 5 decimal places in the OP's linked file, the Lorentzian function evaluation spacing of 0.01 $cm^{-1}$ introduces negligible error.

| improve this answer | |
$\endgroup$
  • $\begingroup$ As usual the OP has vanished despite giving him useful information. $\endgroup$ – M. Farooq Aug 1 at 5:43
  • 1
    $\begingroup$ @M.Farooq Yes. And yet other OPs quickly accept so-so or half-baked answers, before an actual thoughtful answer gets posted. $\endgroup$ – Ed V Aug 1 at 14:26
1
$\begingroup$

The problem (noted in another answer) is that the line broadening function is not centered, therefore the effect of the convolution, in addition to modifying the peak widths, is to shift the spectrum. An off-center Lorentzian (such as used by the OP) is itself a convolution of a centered Lorentzian and a shifted delta function. If you ignore the Lorentzian for a moment, the effect of the shifted delta function is to shift the spectrum.

If a centered LB function is used, as shown in the following figure, the problem is largely resolved:

enter image description here

I did this by using the following convolution function:

L = 1./(1 + 20*(vsim - mean(vsim)).^2);

where vsim are the frequencies in the interpolated spectrum.


Aside

The code in OP leaves much to be desired: not commented, unused blocks of code and variables, non-optimal use of Matlab's intrinsic matrix handling properties and inbuilt functions. However this is not a subject for this SE - perhaps better to ask over at code review or another SE.

The following performs the desired convolution. It is not ideal (see previous comment) but will get the job done:

HCN2nu2 = importdata('./HCN 2nu2.txt');
HCN = HCN2nu2(1:111,1);
AX = HCN2nu2(1:111,2);
v = []; s = []; L = [];
range = 6200:0.01:6700;
p = mean(range);
for i = range
    for j = 1:111
        if i == round(HCN(j),2)
            z = AX(j);
            break
        else
            z = 0;
        end
    end
    v = [v; i];
    s = [s; z];
    t = 1/(1 + 4*(i - p)^2);
    L = [L; t];
end

figure
plot(v,s)
k = conv(L,s, 'same');
hold on
plot(v,k,'k')
| improve this answer | |
$\endgroup$
  • $\begingroup$ BuckThorn, Mathematics never fails to surprise us. I am wondering the why your mean centered Lorentzian also works. Suppose we had one stick at w= 6500, now the convolution property is that the two centroids add. For example, your Lorentzian is 6450 and the stick is at 6500, the resulting convolution is 6450+6500. Not sure, why you do not see a shift. $\endgroup$ – M. Farooq Jul 29 at 14:33
  • $\begingroup$ BuckThorn, What is the source of your 20? Please see the Lorentzian equation (1) which Ed and I am using. $\endgroup$ – M. Farooq Jul 29 at 16:51
  • $\begingroup$ @M.Farooq In essence, by subtracting the mean in the equation shown the Lorentzian is centered in the window, which is exactly what you suggest the OP do to fix the problem. I happened to implement it without seeing your answer first. After I posted my answer I read yours more carefully and realized that other than providing guidance with the code my answer was not going to enlighten anyone on the finer points of FT and convolution, which you already covered. I am tempted to post more code but think the OP has to conduct a substantial code review. $\endgroup$ – Buck Thorn Jul 29 at 19:49
  • $\begingroup$ The 20 is proportional to the inverse linewidth. $\endgroup$ – Buck Thorn Jul 29 at 19:49
  • $\begingroup$ PS The detailed values of the frequency scale shown in the plot of the LB function are not meaningful. What's important is registry with the interpolated stick spectrum, as you also suggest. $\endgroup$ – Buck Thorn Jul 29 at 19:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.