考慮由點(x1,y1), (x2,y2), (x3,y3), ... ,(xn,yn)算法以平滑的曲線，同時保持其下的面積恆定
限定的離散曲線定義一個常數SUM = y1+y2+y3+...+yn。假設我們改變了k個y點的數值（增加或減少），使得這些變化點的總和小於或等於常數SUM。
什麼是最好的方式來調整給出以下兩個條件以外，其它Ÿ點：
SUM。一個簡單的解決方案將是，以限定某些增量如下：
delta = (ym1' + ym2' + ym3' + ... + ymk') - (ym1 + ym2 + ym3 + ... + ymk')
，並在點的其餘部分均勻分佈這個delta。這裏ym1'是修改後的修改點的值，ym1是修改前的修改點的值，以給出修改的總差異。
但是，這並不能確保完全平滑的曲線，因爲靠近變化點的區域會顯得粗糙。這個問題是否存在更好的解決方案/算法？
我已經使用了下面的方法，雖然它有點OTT。
考慮將d [i]添加到y [i]中，以獲得平滑值s [i]。 我們尋求最小化
S = Sum{ 1<=i<N-1 | sqr(s[i+1]-2*s[i]+s[i-1] } + f*Sum{ 0<=i<N | sqr(d[i])}
中的第一項是曲線的（近似）的二階導數的平方和，而第二項懲罰移動從原來的路程。 f是（正）常數。小代數重鑄此作爲
S = sqr(||A*d - b||)
其中矩陣A有一個很好的結構，並且實際上A「* A是五對角線，這意味着該正規方程（即d = INV（A」 * A） * A'* b）可以有效解決。請注意，d是直接計算的，不需要初始化它。
鑑於溶液d到這個問題，我們可以計算出溶液d ^到相同的問題，但與約束在其中'* d = 0（其中一個是所有那些的載體）這樣
d^ = d - (One'*d/Q) * e
e = Inv(A'*A)*One
Q = One'*e
f有什麼價值？一個簡單的方法就是在各種fs的樣本曲線上嘗試這個過程，並選擇一個看起來不錯的值。另一種方法是選擇平滑度的估計值，例如二階導數的均方根值，然後是應該達到的值，然後搜索給出該值的f。作爲一般規則，f越大，平滑曲線將會越平滑。
所有這些的一些動機。目標是找到給定的'平滑'曲線'接近'。爲此，我們需要一個平滑度量度（S中的第一項）和一個接近度量度（第二項，爲什麼這些度量值？每個都很容易計算，並且在變量中每個都是二次的（d []） ;這意味着問題成爲線性最小二乘法的一個實例，其中有高效的算法可用，而且每個和中的每一項取決於變量的附近值，這又意味着'逆協方差'（A' * A）將有一個帶狀結構，所以最小二乘問題可以有效解決。爲什麼引入f？好吧，如果我們沒有f（或將它設置爲0），我們可以通過設置d [i] = -y [i]，獲得完美平滑的曲線s [] = 0，這與y曲線無關。另一方面，如果f是巨大的，那麼爲了最小化s，我們應該專注於第二項，並設置d [i] = 0，我們的「平滑」曲線就是原始的，所以我們有理由假設隨着f的變化，相應的解會發生變化之間非常光滑，但遠離y（小f），接近y但有點粗糙（大f）。
人們常說，我在這裏主張的正則方程是解決最小二乘問題的一種不好的方法，而且這通常是正確的。然而，對於'好'的帶狀系統 - 就像這裏一樣 - 通過使用正態方程的穩定性的損失並不是很大，而速度的增益是如此之大。我已經使用這種方法在合理的時間內使數千個點的曲線平滑。
要了解A是什麼，請考慮我們有4點的情況。然後，我們的表情歸結爲：
sqr(s[2] - 2*s[1] + s[0]) + sqr(s[3] - 2*s[2] + s[1]) + f*(d[0]*d[0] + .. + d[3]*d[3]).
如果我們替代S [i] = Y [I] + d [I]在這方面我們得到的，例如，
s[2] - 2*s[1] + s[0] = d[2]-2*d[1]+d[0] + y[2]-2*y[1]+y[0]
等我們看到，這是SQR（|| A *分貝||），我們應該採取
A = (1 -2 1 0)
(0 1 -2 1)
(f 0 0 0)
(0 f 0 0)
(0 0 f 0)
(0 0 0 f)
and
b = (-(y[2]-2*y[1]+y[0]))
(-(y[3]-2*y[2]+y[1]))
(0)
(0)
(0)
(0)
在實施中，不過，你可能不希望形成A和b，因爲他們只打算用於形成正常方程項，A'* A和A'* b。直接累積這些會更簡單。
怎麼樣保持相同的動態範圍？
min0,max0y - 值min1,max1y - 值線性插值的所有值以匹配原始最小值最大值y
y=min1+(y-min1)*(max0-min0)/(max1-min1)
就是它
不知道的領域，但這應該保持形狀更接近原來的。我得到了這個想法，現在在閱讀你的問題，我現在面對類似的問題，所以我儘量代碼，並馬上想試試+1的讓我這個想法:)
可以適應這一點，並與區域相結合
因此，在此之前，計算該區域並應用＃1 ..＃4並在那之後計算新的區域。然後將所有值乘以old_area/new_area比率。如果你也有負面的價值，而不是計算絕對面積，那麼你必須分別處理正面和負面的區域，並找到最適合攤位的原始面積的倍增率。
[EDIT1]一些結果恆定動態範圍
正如你可以看到形狀稍微移動到左邊。每張圖片在應用了幾百次平滑操作之後。我想到的是細分到地方最小值最大值間隔，以提高這一點...
[EDIT2]已經完成了過濾器用於我自己的目的
void advanced_smooth(double *p,int n)
{
int i,j,i0,i1;
double a0,a1,b0,b1,dp,w0,w1;
double *p0,*p1,*w; int *q;
if (n<3) return;
p0=new double[n<<2]; if (p0==NULL) return;
p1=p0+n;
w =p1+n;
q =(int*)((double*)(w+n));
// compute original min,max
for (a0=p[0],i=0;i<n;i++) if (a0>p[i]) a0=p[i];
for (a1=p[0],i=0;i<n;i++) if (a1<p[i]) a1=p[i];
for (i=0;i<n;i++) p0[i]=p[i]; // store original values for range restoration
// compute local min max positions to p1[]
dp=0.01*(a1-a0); // min delta treshold
// compute first derivation
p1[0]=0.0; for (i=1;i<n;i++) p1[i]=p[i]-p[i-1];
for (i=1;i<n-1;i++) // eliminate glitches
if (p1[i]*p1[i-1]<0.0)
if (p1[i]*p1[i+1]<0.0)
if (fabs(p1[i])<=dp)
p1[i]=0.5*(p1[i-1]+p1[i+1]);
for (i0=1;i0;) // remove zeros from derivation
for (i0=0,i=0;i<n;i++)
if (fabs(p1[i])<dp)
{
if ((i> 0)&&(fabs(p1[i-1])>=dp)) { i0=1; p1[i]=p1[i-1]; }
else if ((i<n-1)&&(fabs(p1[i+1])>=dp)) { i0=1; p1[i]=p1[i+1]; }
}
// find local min,max to q[]
q[n-2]=0; q[n-1]=0; for (i=1;i<n-1;i++) if (p1[i]*p1[i-1]<0.0) q[i-1]=1; else q[i-1]=0;
for (i=0;i<n;i++) // set sign as +max,-min
if ((q[i])&&(p1[i]<-dp)) q[i]=-q[i]; // this shifts smooth curve to the left !!!
// compute weights
for (i0=0,i1=1;i1<n;i0=i1,i1++) // loop through all local min,max intervals
{
for (;(!q[i1])&&(i1<n-1);i1++); // <i0,i1>
b0=0.5*(p[i0]+p[i1]);
b1=fabs(p[i1]-p[i0]);
if (b1>=1e-6)
for (b1=0.35/b1,i=i0;i<=i1;i++) // compute weights bigger near local min max
w[i]=0.8+(fabs(p[i]-b0)*b1);
}
// smooth few times
for (j=0;j<5;j++)
{
for (i=0;i<n ;i++) p1[i]=p[i]; // store data to avoid shifting by using half filtered data
for (i=1;i<n-1;i++) // FIR smooth filter
{
w0=w[i];
w1=(1.0-w0)*0.5;
p[i]=(w1*p1[i-1])+(w0*p1[i])+(w1*p1[i+1]);
}
for (i=1;i<n-1;i++) // avoid local min,max shifting too much
{
if (q[i]>0) // local max
{
if (p[i]<p[i-1]) p[i]=p[i-1]; // can not be lower then neigbours
if (p[i]<p[i+1]) p[i]=p[i+1];
}
if (q[i]<0) // local min
{
if (p[i]>p[i-1]) p[i]=p[i-1]; // can not be higher then neigbours
if (p[i]>p[i+1]) p[i]=p[i+1];
}
}
}
for (i0=0,i1=1;i1<n;i0=i1,i1++) // loop through all local min,max intervals
{
for (;(!q[i1])&&(i1<n-1);i1++); // <i0,i1>
// restore original local min,max
a0=p0[i0]; b0=p[i0];
a1=p0[i1]; b1=p[i1];
if (a0>a1)
{
dp=a0; a0=a1; a1=dp;
dp=b0; b0=b1; b1=dp;
}
b1-=b0;
if (b1>=1e-6)
for (dp=(a1-a0)/b1,i=i0;i<=i1;i++)
p[i]=a0+((p[i]-b0)*dp);
}
delete[] p0;
}
所以p[n]是輸入/輸出數據。有跡象表明，可以調整等幾件事情：
<0.8,0.8+0.35/2>）（在for循環現在5）1.0意味着沒有變化的主要思想是背後：
計算權重，用於平滑
所以鄰近局部極值是輸出
光滑
[注意事項]
我也嘗試恢復該地區，但是這是與我的任務不兼容，因爲它扭曲的形狀很多。所以，如果你真的需要這個區域，那麼就關注那個區域，而不是那個形狀。平滑造成信號大部分縮小後區域恢復的幅度上升。
實際過濾器狀態沒有可標記的形狀側移（這是我的主要目標）。更多的顛簸信號某些圖像（原過濾器是在這個可憐的極端）：
正如你可以看到沒有明顯的信號形狀轉移。當地有極端傾向很重的平滑之後形成鋒利的尖刺但預計
希望它可以幫助...
我喜歡這個。它具有最小能量樣條的概念。 Upvoted它。 – fang 2014-12-04 19:09:31
我不明白解決方案。我們爲什麼要最小化S和A？其實我知道發生了什麼，但是爲什麼它似乎完全讓我昏迷。 – Sohaib 2014-12-05 05:22:20
很好地解釋了很多。只是更清楚一點，那麼我會將其標記爲正確。d的初始值應該是多少？ 又是什麼？你能否解釋一下你是如何得到第二個方程的：S = sqr（|| A * d-b ||）'。 – Sohaib 2014-12-06 16:02:48