my Butterworth lowpass formulas do not agree with Fisher webpage I want to implement Butterworth low-pass filter. Thanks to this question, I have found out that the filter coefficients can be generated using Tony Fisher web-site or using his code. But the problem arose when I had tried to verify his formulas myself.

Wikipedia says that the derivation of low-pass formulas is simple: we start with the Butterworth polynomial of order n. I took n=1
$$B_1\left(\frac{s}{\omega_{cutoff}}\right)=1+\frac{s}{\omega_{cutoff}}$$

(note that $\omega_{cutoff}$ is angular frequency, not usual one) then do bilinear transform
$$s = 2f_{sampling}\cdot\frac{z-1}{z+1}$$

($f_{sampling}$ is usual frequency, not angular) and rewrite the resulting fraction in the form with non-positive powers of $z$.

To make the story short, my final formula for transfer function is
$$H(z) = \frac{Y(z)}{X(z)}$$
where
$$Y(z) = \frac{a}{1+a}+\frac{a}{1+a}z^{-1}$$
$$X(z) = 1-\frac{1-a}{1+a}z^{-1}$$
and
$$a = \frac{\omega_{cutoff}}{2f_{sampling}}$$
and the resulting formula is
$$y[n] = \frac{a}{1+a}x[n]+\frac{a}{1+a}x[n-1]+\frac{1-a}{1+a}y[n-1]$$
For testing I use $f_{cutoff}=1$ (which is $\omega_{cutoff}=2\pi$) and $f_{sampling}=30$.

We should not worry about coefficients in front of $x[]$ since Fisher multiplies them by gain factor anyway, but the coefficient in front of $y[n-1]$ equals

$$\frac{1-a}{1+a} = 0.8104139027$$

while Fisher's web-page for $f_{cutoff}=1$ and $f_{sampling}=30$ gives $0.8097840332$.

If you had a patience to finish this reading, may be you can explain, where I (or Fisher?) am wrong.

The problem is in the way you apply the bilinear transform. You have to use the appropriate (pre-)warping of the frequencies. Since the bilinear transform warps the frequency axis, you have to make sure that the corner frequency of the discrete-time filter is correct. One way to do that is as follows. The bilinear transform is defined as

$$s=k\frac{z-1}{z+1}\tag{1}$$

with some constant $k$ yet to be defined. If we denote the analog frequency by $\Omega$, and the normalized discrete-time frequency by $\omega$, Eq. (1) becomes for $s=j\Omega$ and $z=e^{j\omega}$

$$j\Omega=k\frac{e^{j\omega}-1}{e^{j\omega}+1}=k\frac{e^{j\omega/2}-e^{-j\omega/2}}{e^{j\omega/2}+e^{-j\omega/2}}=jk\tan(\omega/2)\tag{2}$$

Eq. (2) describes the frequency warping caused by the bilinear transform. If we use an analog lowpass filter with a normalized corner frequency $\Omega_0=1$, we must choose the constant $k$ such that for the desired discrete-time corner frequency $\omega_0$ the term on the right-hand side of Eq. (2) becomes $1$:

$$k=\frac{1}{\tan(\omega_0/2)}\tag{3}$$

where $\omega_0$ is the desired corner frequency of the digital filter. Eq. (3) and Eq. (1) define the appropriately normalized bilinear transform that you must use.

So for your example, the normalized first-order analog Butterworth lowpass transfer function is given by

$$H(s)=\frac{1}{1+s}\tag{4}$$

Applying the bilinear transform gives

$$H(z)=\frac{1}{1+k\frac{z-1}{z+1}}=\frac{z+1}{z(1+k)+1-k}=\frac{1}{1+k}\cdot\frac{1+z^{-1}}{1+\frac{1-k}{1+k}z^{-1}}\tag{5}$$

Ignoring the gain term $1/(1+k)$, the corresponding difference equation is

$$y[n]=x[n]+x[n-1]-\frac{1-k}{1+k}y[n-1]\tag{6}$$

With a desired corner frequency $\omega_0=2\pi/30$ you get from (3) $k=9.5144$ and $(1-k)/(1+k)=-0.80978$, just like on Fisher's website.


How do you compute for fc and t in modulating an I and Q component I know that in a $2^b$-ary QAM system you need to modulate the $I$ and $Q$ components through the mixer and are separated by 90° through a local oscillator (by applying a $\cos$ and $\sin$).

By mixing them through the mixer, the $I$ component will be multiplied by a $\cos(2\pi f_c t)$ and $Q$ by $\sin(2\pi f_c t)$. My question is, if I have a 8 bit sampler and the sampling frequency I used is 10 kHz, what will be my value for the carrier frequency $f_c$ and the time sweep, $t$?

My idea is that I can compute the carrier frequency by computing for the symbol rate (bit rate over the number of symbols) and multiplying it by 2 to satisfy the Nyquist Sampling rate. Is that true? And I really don't have any idea what increment value of $t$ I will use.

The answer is: it depends on the pulse shape.

A 16-QAM baseband signal can be written as $$s_{BB}=\sum_{k=0}^{N-1} a_k p(t-kT_p),$$ where $a_k$ is a complex symbol taken from your 16-QAM constellation, $R_p=1/T_p$ is the pulse rate (how many pulses per second you transmit), and $p(t)$ is the pulse shape (ideally, $p(t)$ is a Nyquist pulse).

The bandwidth of $s_{BB}(t)$ is the same as the bandwidth of $p(t)$ (assuming you're transmitting uncorrelated data). Assume this bandwidth is $B$. Then, after mixing with the carrier, the passband signal's maximum frequency is $f_c+B$, and you need a sampling frequency $f_s>2(f_c+B)$ to satisfy Nyquist.

If you don't know what the pulse shape is, then you're probably transmitting rectangular, non-overlapping pulses. As a rule of thumb, I define the bandwidth of a rectangular pulse as five times the pulse rate: $B=5R_p$. Note that many textbook authors, such as Stallings, use $B=R_p$. I prefer to be conservative and define a larger bandwidth.

Then, your sampling frequency is $f_s>2(f_c+5R_p)$.

As an example, let's say that you want to transmit 1000 pulses per second, with a carrier frequency of 1 MHz, and rectangular pulses. Then,

fs = 2*(1e6+5*1000);  % or larger
t = 0:1/fs:Tend       % Tend depends on how many pulses you transmit


Spectral Shape and properties I started learning audio dsp and the book I am working through is loaded with formulas and most of the time expects the reader to know what the relation is between these functions and why they are being used in the first place.

I'd reading a section on spectral shape, which includes sub sections on rolloff, flux, centroid, spread etc.

On what data are the functions actually applied starting from the audio source?

Im trying to get a grip on the actual process behind extracting information from audio.

The book also mentions

  "It is also common to use the power spectrum instead of the magnitude
  spectrum"

How does one change between spectrums and when is it done in the extraction process?

Are they all applied to data in the frequency bins after running a Fourier transform?

For spectral features like roll-off, flux, centroid, etc, the process is as follows:

Before delving into the extraction process, a pre-processing phase is sometimes necessary, this includes:
* Converting stereo signals to mono by summing both channels.
* Resampling the signal to a target sample rate (if a collection of signal of heterogenous sample rates are used).

Now, we can assume that we are dealing with a mono signal of known sample rate.

At each time $t = 0, \Delta, 2\Delta...$ at which the feature has to be calculated (for example, $\Delta = 0.04s$ if 25 values of the feature are to be extracted per second of audio):

Extract the block of signal in the range $[t, t+l]$. $l$ is the analysis window duration and must be greater or equal than $\Delta$. There is a compromise to be found between using short $\Delta$ and $l$ (good temporal resolution) and using a long $l$ (good frequency resolution).
Multiply this signal chunk by a windowing function such as the Hann window. This prevents spectral leakage which might cause unwanted results (for example, report abnormally high values of spectral spread for a purely sinusoidal signal).
Eventually, pad the resulting signal with zeros. The purpose is twofold: make sure the number of samples processed at the next step is a power of 2 (in case the FFT library you use is rustic and does not handle other sizes); and in some applications, artificially boost the frequency resolution through spectral interpolation.
Compute the FFT of the windowed, padded, block of signal.
Discard the upper half of the coefficients returned by the FFT function (they correspond to "negative frequencies" and are of no interest).
Compute either the magnitude or square magnitude of the remaining coefficients. This gives the spectrum on which spectral features are computed.
From there, use whatever formula describes the feature to be extracted. The formulas usually involve $A(k)$ (squared magnitude or magnitude in the $k^{th}$ FFT bin) and $f(k) = \frac{k . \text{sample rate}}{\text{FFT size}}$ (Center frequency of the $k^{th}$ FFT bin).

There is little consensus in the audio signal processing community as to whether the magnitude or power spectrum should be used for spectral features - I have routinely seen both used in research code.


Intuition behind lattice filters In papers and textbooks and on Wikipedia, I regularly stumble upon so-called "lattice filters". They always look horribly more complex than their 'non-lattice counterpart' and I usually don't really understand what's going on (and why). I could not really find a useful introduction to lattice filters and the fundamental idea behind them. What I believe to have understood (please correct me if I'm wrong) is that they represent some form of very efficient implementation scheme for classical filters, e.g., RLS?
Can somebody explain the intuition behind lattice filters? What are they, which problem do they solve, and how do they solve it? References to readable introductions are also highly appreciated.
(Note that I have actually had a fair amount of DSP in my undergrad studies, but lattice filters were never even mentioned.)

The intuition behind lattice filters probably is connected with the propagation of waves and modeling wave propagation with discrete cells.


Transfer function estimation using system identification I am posting a control system problem here because of this post reference.
Question:
I have input and output datasets. These are graphically shown below.

The blue line is output and the pink line is input which was applied during specific time intervals. This dataset has been obtained from a machine that has an existing manually tuned the PID controller.
I want to make an LQI or more perfectly tuned PID controller using bode plots. However, to do that, I will have to identify the transfer function or corresponding state space model.
What is the detailed procedure to obtain the transfer function from the above data set (because this data has been obtained from a machine which already has a manually tuned PID controller controlling it)? The manually tuned values are $K_p=2.2, K_i=0.01, K_d=2$.

My attempt
I used a system identification application and obtained a state space model with a $75 \%$ fit. I tried implementing a PID controller for it and it required large $K_p, K_d, K_i$ values for the required design parameters.

The gain values used were

These cannot be implemented on a system since it has comparatively very low $K_p, K_i, K_d$ values (manually tuned).

Let us denote the output feedback SISO plant transfer function $\frac{Y(s)}{U(s)}=P(s)$ so that the open-loop (OL) transfer function of the system loop is $CP$ where $C(s)=K_p + K_d s + \frac{K_i}{s}$ and the closed-loop (CL) transfer function is $\frac{CP}{1+CP}$. From the attempt at the solution to the problem in the OP, we can see that a second order model might be considered a good idea.
Since the controller transfer function ($K_p, K_d s, {K_i}$) and the system is using direct output feedback, the control input (or actions) may be obtained from the provided data as $u(t) = K_p y(t) + K_d \dot{y}(t) + K_i \int_0^t y(t) \; dt$. The data available is then $\{y(t_i), u(t_i)\}$ which can be directly fed into a system identification tool. Alternately one may used least squares to obtain the second order model parameters ($A, B)$ to fit the state transition model ${\Huge[} \begin{matrix} y \\ \dot{y}\end{matrix}{\Huge]} = A {\Huge[} \begin{matrix} y \\ \dot{y}\end{matrix}{\Huge]} + B {\Huge[} \begin{matrix} 0 \\ u\end{matrix}{\Huge]}$. Care must be taken not to use the CL data to fit the model instead since, as clearly seen from the model, there is no integral term in the state space representation the model is attempting to use. Finally, from a preliminary observation of the graphical representation of the data, the data can be divided into three distinct data set split by time-intervals in which the system undergoes large control actuations which seem to become saturated. These three data sets are observed to have the typical second order behavior reflecting a first or second order system controlled by a PID controller.


Object-Recognition using SIFT: Switching the order of images creates two different results I've developed a large test-system which can use every interest-point-detector and -descriptor. I use it for evaluation of the given algorithms with the standford-image-collection of cd-covers. This post will not go any deeper with any source-code. I only found out that switching to images (no referred as firstImg and secondImg) will create to different results. The Ransac-Threshold is used to compute the homography and epsilon is the reprojection-error. For this example I use the Scale-Invariant-Feature-Transform by David Lowe.

First test with the secondImg as the queue-image creates the following results I've written down:

Inliner: "1672"
Number_of_Keypoints_FST: "368"
Number_of_Keypoints_SND: "1749"
Outlier: "77"
Ransac: "65"
Epsilon: "10"

Now I'm switching the order. The firstImg is the queue.

Inliner: "310"
Number_of_Keypoints_FST: "1749"
Number_of_Keypoints_SND: "368"
Outlier: "58"
Ransac: "65"
Epsilon: "10"

Questions: I'm confused that switching the order of the images affects the results. Should the matched keypoints be identical for the matched ones when switching the order of the images? Should the inliner be identical?

ARGH. I feel ashamed. How about a cross validation? That's the solution. Now there are almost exact matches ... almost. =) Thank you, Libor. The second comment gave me a hint. =)


Correction of amplitude after zero padding for upsampling purposes I have time sequence in which the data is sampled at 0.8 Hz. The data is related to chromatography (chemical analysis), that is why the sampling frequency is relatively low. The instrument cannot sample faster at this moment.
I was exploring the idea of upsampling by zero padding in the frequency domain as follows in MATLAB.
FFT_S=fft(Signal); % FFT of Signal of Interest

FFT_ZP=[FFT_S(1:length(t)/2,1); zeros(1000,1); FFT_S(length(t)/2+1:end,1)]; % Zero padding with 1000 zeros.

Signal_Up=real(ifft(FFT_ZP)); % Upsampled data

The original data consists of 716 points. The upsampled data has 1716 points but the amplitude has reduced - which is undesirable.
Is there a simple multiplying factor to correct the amplitude in MATLAB based on total number of points before and after upsampling?
EDIT:
Fortunately for this analytical chemistry purpose the trade offs of zero padding are not relevant. Qualitatively, in order to keep the amplitude the same, I found that if we double the sampling rate, the amplitude after inverse fft had to be multiplied by 2, if we triple the sampling rate by zero padding, the amplitude had to be multiplied by 3. Ignoring the trade offs, there must be a generalized method to correct the final amplitude based on the initial and final number of data points?
Thanks.

If you use an 1/N normalized DFT, then you shouldn't have any problems with your amplitude when you take the inverse DFT (with no normalization factor/factor of 1).  Consider the case of a pure tone with a whole number of cycles in the frame.  Only one bin pair will be non-zero.  No matter how many extra zeroes you insert, the inverse DFT will reconstruct the same signal in the same duration, just appropriately upsampled.  The amplitude of the signal will be twice the magnitude of of each bin value independent of your sample count.

If the normalization factors are already hard coded in your routines as "1" forward, "1/N" inverse, you can simply use them in reverse roles.  Take the "inverse" of your signal, zero pad it (splitting the Nyquist), then take the "forward".
Alternatively, just multiply your final results by a factor of $\frac{N_{new}}{N_{old}}$, which confirms your observations in the update.


Contradiction between complex baseband and real-valued baseband Suppose a direct-conversion transmitter+receiver (ideal transmitter) along with its filters. In complex baseband, the filters in the signal chain can be modeled as a complex-valued FIR filter ($y$, $c$, $x$ are complex-valued!):

$$
y = \sum_{k=0}^K c_k x[n-k]
$$

Rewriting this in cartesian coordinates gives:

$$
y_i + j y_q = \underbrace{\sum_{k=0}^K \left(a_k x_i[n-k] - b_k x_q[n-k]\right)}_{y_i} + j \sum_{k=0}^K \left( b_k x_i[n-k] + a_k x_q[n-k] \right)
$$

Now only considering the real part, it can be seen that it consists of two different filters with which the real and imaginary input signals are filtered.

Now consider the same system from a practical perspective: After the DAC and upconversion, the RF signal looks as follows:

$$
x_{rf} = x_i \cos \omega_c t + x_q \sin \omega_c t
$$

The output sequence $y_i$ is obtained by multiplying $x_{rf}$ with $\cos(\omega_c t +\phi)$ and filtering it with a filter called $H$:

$$
y_i = H( x_{rf} \cos(\omega_c t +\phi) ) \\
\approx H(x_i/2 \cos\phi + x_q/2 \sin\phi) \\
= \frac{\cos\phi}{2} H(x_i + x_q \tan\phi)
= G(x_i + x_q \tan\phi)
$$

The filter $H$ (or $G$) is modeled as a FIR filter:

$$
y_i = \sum_{k=0}^K \left( g_k x_i[n-k] + g_k \tan\phi x_q [n-k] \right)
$$

From the equation above it can be seen that the real and imaginary part are filtered through only one filter (they differ just by the constant $\tan\phi$ !).

Where does this contradiction come from?

PS: I know that the first approach is the correct one because it gives the correct results. I do not understand why it is not consistent with my second approach

As a sidenote:

When you demodulate your signal $x_{\mathrm{rf}}(t)$  you say that there is some phase difference $\phi$. This translates to crosstalk (i.e., the quadrature and in-phase components influence each other at the receiver). It's easy to see that if $\phi = 0$ (namely, we know exactly when to sample the received signal), there is no crosstalk at all:

$$y_i = \sum_{k=0}^K \left( g_k x_i[n-k] + g_k \tan(0) x_q [n-k] \right) =  \sum_{k=0}^K  g_k x_i[n-k]$$

In the baseband domain, this means that the filter $c_k$ is purely real (if it was purely imaginary, then there wouldn't be crosstalk either, you just have to invert the channels to retrieve the information correctly).

The answer:

Whether there is crosstalk or not, both of your equations get to the same result, but with different names. Using your notation, $a_k$ in the first equation would be equivalent to $g_k$ in the second one, and $b_k$ to $-g_k\tan(\phi)$. I believe that your confusion arose from thinking of $a_k$ and $b_k$ as two different filters (which is correct), instead of thinking of them as the real and imaginary part of a single filter ($c_k$).


Recognizing object of interest in object tracking using Lucas-Kanade method I am able to find out the velocities of pixels using the Lucas-Kanade method. Now I want to use it for tracking. I am using a Kalman filter taking $x$, $y$, $v_x$ and $v_y$ ($(x,y)$ is position of a particular object pixel and $(v_x,v_y)$ is velocity of a pixel) as state variables of the Kalman filter. Now the problem I am facing is how will I get to select the pixels that are representing the object to be tracked?

Tracking objects can be simple to complex depending on what type of background you have, whether the background is static or moving, whether the object is clearly distinguishable or can share similar properties with background there by blending in the background, and whether we have unique object at hand or we are tracking multiple objects which may cross the trajectories.

The good news is, the solution to all of these exists. Depending on the nature of task at hand, you may choose simple to complex solution.

What exactly we need?

As you implement any tracking algorithm, you have an estimate of where the new position of the object is likely be. Now you need to verify that indeed the object is at the same position or different. The difference is the error that you feedback to your model. For this you need to identify/segment the object in the vicinity where your projected value. In other words you need to detect that object with in a given window.

If the object's appearance or key behaviour doesn't change quite drastically you can capture feature vectors from original reference and compare the same around the vicinity where your target is expected to be. This is search of feature vector within a window based on some distance criteria.

So the question boils down to - what should be the feature vector and what should be the distance matrix?

There are many possible answers. Writing down some of them (from simple to complex) with remark on when they might work.

Color based features: This include identifying mean or dominant color, color histogram over a local window etc. If the object color is unique and doesn't blend with background then you can identify object with this manner. See [Ref 1] and [Ref 2] below
Texture based features: There is a whole range of features here. The famous Gabor filter banks, then there is co-occurrence matrix, you can take gradient of images or find edges and have edge histograms over local regions, or you can identify Local binary patterns near the regions.
See [Ref 3]. Also, many approaches combines Color and Texture information for more effectiveness. See [Ref 4] and [Ref 5]
Shape based features: Shape based technique works when the object you are tracking has a particular shapes that distinguishes from other objects and from background. This includes many techniques. See [Ref 6], [Ref 7] and [Ref 8] below.

More complex features:

If the above features doesn't work out, specifically if the object is complex and the background as well -then we can use these features:

SIFT - The scale-invariant feature transform allows to compute key points which characterises the object and hence you can compare the same object from past frame against a new frame.  See [Ref 9]
SURF, FAST, ORB:  - If SIFT turns out to be blocked because of patent issue, and due to it's speed, then you can use it's alternatives which are SURF, FAST, ORB. See [Ref 10] for a good comparison.
HoG and Haar: Finally you can use HoG - which is a good generic object detection framework. [Ref 11] and Haar based classifiers when the problems are really complex. [Ref 12]

I am not prescribing any specific solution, and none of the above could be your direct solution - yet this will help you thinking through the problem at hand. You can pick whichever feature works for you irrespective of tracking methods that are used in those respective papers.

Go through the paper Good Features to Track  thoroughly and it will give you much better perspective on the same.

What about distance criteria?

For each different feature typically different distance measure should make sense. But, distance matrix is usually dependent on your feature vector itself.

[Ref 1]: Object Tracking Based on Color Information
Employing Particle Filter Algorithm Budi Sugandi, Hyoungseop Kim et. al.

[Ref 2]: Applying a New Spatial Color Histogram in Mean-Shift
Based Tracking Algorithm Dong Xu, Yimin Wang, and Jinwen An

[Ref 3]: Texture Features for Browsing and Retrieval of Image Data B.S. Manjunathi and W.Y. Ma

[Ref 4]: ROBUST OBJECT TRACKING USING JOINT COLOR-TEXTURE HISTOGRAM JIFENG NING, LEI ZHANG and DAVID ZHANG, CHENGKE WU

[Ref 5]: Multi-Object Tracking Using Color, Texture and Motion Valtteri Takala and Matti Pietikainen

[Ref 6] :  A Survey on Feature Extraction Techniques for Shape based Object Recognition Mitisha Patel, Purvi Tandel

[Ref 7]: RAPID AND ROBUST HUMAN DETECTION AND TRACKING BASED ON OMEGA-SHAPE
FEATURES Min Li, Zhaoxiang Zhang, Kaiqi Huang and Tieniu Tan

[Ref 8]: An Evaluation of Local Shape-Based Features for Pedestrian Detection Edgar Seemann, Bastian Leibe, Krystian Mikolajczyk, Bernt Schiele

[Ref 9]: Object Tracking using SIFT Features and Mean Shift Huiyu Zhou, Yuan Yuan, and Chunmei Shi

[Ref 10]: Image Matching Using SIFT, SURF, BRIEF and ORB: Performance Comparison for Distorted Images Ebrahim Karami, Siva Prasad, and Mohamed Shehata

[Ref 11]: Part-based Multiple-Person Tracking with Partial Occlusion Handling Guang Shu, Afshin Dehghan et. al.

[Ref 12]: Real-time hand tracking using a set of cooperative classifiers based on Haar-like features Andre L. C. Barczak and Farhad Dadgostar


Why MATLAB only allows cpsd and mscohere between vectors with the same length? I am studying the applications of finding coherence in the frequency spectrum between two signals for pattern recognition. I have a recording that is way longer than the sampled signal that I want to compare. However, MATLAB only allows cpsd and mscohere between vectors with the same length. I figured that there must be a good reason for that, but I could not find it.

Also, I was digging through the math of the procedure and I figure that I could do a crosscorrelation between the two (although it would be raw, i.e., not unbiased), and then apply the Fourier transform to obtain the correlated frequencies. Do you think it is sound or make sense to do such a thing?

The calculations are done in the frequency domain using pairs of FFT's which is why the vectors need to be the same length.  You could truncate the longer vector to the length of the shorter. If the longer is much longer, you could select multiple (different) truncated segments.

You can do things in the time domain like cross correlation, but mscohere and cpsd average (the default of 8 segments).  A single crosscorrelation isn't equivalent.

When reading "does it make sense" (you're not the only one who asks), Internally, I think this is a lack of confidence and I don't want to discourage this person from trying, but this is whacko .  Try it. MATLAB is easy.  It's designed to explore the whacko.


How many samples per cycle(time period) is actually enough to reproduce the sine wave sound properly? I have been working on Active Noise Cancellation in real time for Air vent in kitchen.
Initially for the demonstration purpose, I am using sinusoidal sound ( 40k Hz Fs,1000 samples, 2 kHz frequency, 1 Amplitude). I am using simple LMS Algorithm initially to test this, I could produce the anti noise with the same frequency and amplitude but because it is a real time application, I could be able to produce the sine wave with 30-31 samples per each cycle of sine wave.

The problem is the sound coming through the noise cancelling speakers is not good compared to the source noise, there is some buzz noise coming in addition (kind of background noise). This could mainly be because of distortions in speaker, I don't have any idea how to control this !!
But I am thinking this may be could be because of very less number of samples that I am using for Cancelling speaker.May be the reproduction of the noise is not so good because of less number of samples.

So my question is how many samples per cycle are ideally required for good reproduction of sine wave sound??

It depends on how good your interpolation between samples is.  If you have really good interpolation, anything more than 2 samples per cycle will suffice.  This is not just a theoretical fact, it is my experience in practice using polyphase interpolation that combines 64 adjacent table-lookup samples.

If you're using linear interpolation (which combines 2 adjacent table-lookup samples in the interpolation), I would use no less than 1024 samples per cycle for a high-quality sine wave (about 120 dB S/N).


FFT result of windowed cosine function The frequency spectrum of cosine function is as below.

While adopting a square window function to the cosine function, it is same as implementing convolution of their frequency spectrum as below. The thick dash line is the sum of them. Is it correct to say the points on this dash line are the FFT result of corresponding frequency? I feel it is not true, but don't know where the problem is.

There is no problem in your plot if by fft you mean DTFT (or Discrete Time Fourier transform) and it is true to say "the points on this dash line are the DTFT result of corresponding frequency":
If your function is a cosine:
$$
x[n] = \cos(\omega_0n) \Longrightarrow X(e^{j\omega}) = \pi\delta(\omega-\omega_0)+\pi\delta(\omega+\omega_0)\ ; \qquad|\omega|<\pi
$$
and your window function is defined as follow:
$$
w[n] = \left\{
        \begin{array}{ll}
            1 & \quad |n| \leq N \\
            0 & \quad |n| > N
        \end{array}
    \right.
\Longrightarrow W(e^{j\omega}) = \frac{\sin\left(\omega(\frac{N}{2}+1)\right)}{\sin\left(\frac{\omega}{2}\right)}
$$
as you mentioned if we multiply two functions in time domain, it's equal to periodic convolution in frequency domain:
$$
y[n] = x[n]w[n]\\ \Longrightarrow Y(e^{j\omega})=X(e^{j\omega})\circledast W(e^{j\omega})=\frac{1}{2\pi}\int_{2\pi}^{}X(e^{j\omega})W(e^{j(\omega-\theta)})  d\theta
$$
in this case (in one period):
$$
Y(e^{j\omega})=\frac{1}{2}\frac{\sin\left((\omega-\omega_0)\left(\frac{N}{2}+1\right)\right)}{\sin\left(\frac{\omega-\omega_0}{2}\right)}+\frac{1}{2}\frac{\sin\left((\omega+\omega_0)\left(\frac{N}{2}+1\right)\right)}{\sin\left(\frac{\omega+\omega_0}{2}\right)}
$$
if you plot this function in one period for $\omega_0=\frac{\pi}{5}$ and $N=2$ it look like this:

but if you change the window length to $N=4$ the result changes to:

and the cause for this, is that two main lobes overlap with one another.


Digital delay effect: avoiding clicks on delay time change I'm implementing a digital audio delay effect using a simple circular delay line. Changing the delay time by making the delay line longer or shorter introduces discontinuities in the audio signal, resulting in clicks.

What would be the best (most pleasing for the ears) strategy to avoid such clicks when changing delay time?

If you are not changing the delay length very often, and you don't want to have a Doppler effect that comes from continuously changing the delay length, then try a cross-fade. Both delay lengths should be running simultaneously for a moment, and you would fade the old one out while fading the new one in.


What causes this type of distortion in a constellation I/Q plot? I'm transmitting QPSK through a digital pulse shaping filter (RRC, 35% excess BW, 2 samples per symbol, 64 taps) at a rate of 56Msym/sec.
Receiving it on standard lab equipment (Vector signal transceiver), where it is demodulated with RRC MF and plotting the constellation. At high symbol rates, the phenomenon below starts appearing in the constellation where the points get essentially binned into 4 different locations (looks like some exotic modulation you might call qpsk^2)
One thing of note, this effect is easily taken care of by the equalizer on the RX and reduces the constellation to a much more ideal, low EVM one. Given this, I suspect the effect is related to ISI. One possible thought I had was for example, on the I channel each symbol is perturbed by a neighbor (which could be a + or - symbol) so the effect is two offset bins. Same for Q channel, resulting in the 4 shown. Not sure why this contribution would only come from a single symbol as opposed to many though - in which case I would expect more of a cloud/awgn like effect from adding up contributions of many symbols)
For reference the basic D/A chain looks like:
Generate QPSK IQ -> Filter with RRC -> Halfband Upsamplers -> DAC -> Analog Reconstruction

I don’t think this is necessarily additionally introduced ISI (beyond the ISI of the pulse shaping filter itself, which is zero when there is no timing offset), but may be the result of timing offset error. All the points within a small offset from the ideal sampling location appear to be selected, and the overall pattern of the full constellation showing all samples does not appear to show any distortion.
My theory is we are seeing the results of the overshoot and undershoot you would see in proximity to the ideal sampling location due to the expected response of the pulse shaping filter. What I am describing is better explained with reference to the eye diagrams on the I and Q channel, in particular the case I give below with the higher alpha value where it is quite clear how 4 clusters could result with a small time offset. Notice if the sample is selected at Time = 0 that all the samples will be at the same four I, Q locations, but if we shift the sample time slightly higher there will be 16 possible I Q sample combinations (4 on I and 4 on Q):

Here is the same thing with alpha = 0.25 and the clusters would migrate to a larger cloud:

The OP's case is using alpha =  0.35, so jury is out if this is indeed due to timing offset alone. I recommend reviewing the data with eye diagrams as I have done here to further confirm root cause between residual ISI at the correct sampling location, or timing offset error, or both.


Savitzky–Golay filter vs. IIR or FIR linear filter A traditional IIR / FIR filter (lowpass to remove the high freq oscillations), e.g. moving average,
or a Savitzky-Golay filter

can all be useful to smoothen a signal, such as an envelope signal:

For which application would a Savitzky-Golay filter be more interesting than a classical lowpass?

What makes it different than a standard filter, and what does it add in comparison to standard filters?

Does it adapt itself to the input data?

Is it better for transient preserving?

Have you ever been in an engineering situation one day when you decided "Let's use a SG filter instead of moving average or another FIR lowpass! It's better because this and this and this..."? Then this question is for you!

Since the discussion in the existing answers and comments has mainly focused on what Savitzky-Golay filters actually are (which was very useful), I will try to add to the existing answers by providing some information on how to actually choose a smoothing filter, which is, to my understanding, what the question is actually about.

First of all, I'd like to repeat what has become clear in the discussion spawning from the other answers: the categorization of smoothing filters in the question into linear and time-invariant (LTI) FIR/IIR filters on the one hand, and Savitzky-Golay filters on the other hand is misleading. A Savitkzy-Golay filter is just a standard FIR filter designed according to a specific criterion (local polynomial approximation). So all filters mentioned in the question are LTI filters.

The remaining question is how to choose a smoothing filter. If computational complexity and/or memory are an issue, IIR filters may be preferable over FIR filters, because they typically will achieve comparable noise suppression (i.e., stopband attenuation) with a much lower filter order than FIR filters. But note that if real-time processing is necessary, one possible disadvantage of IIR filters is that they cannot have an exactly linear phase response. So the desired signal will suffer some phase distortions. For offline-processing, the phase distortions can be avoided, even with IIR filters, by applying zero-phase filtering.

Apart from the considerations discussed in the previous paragraph, it is mainly the design criterion that matters, not so much if the filter is FIR or IIR, because any (stable) IIR filter can be approximated with arbitrary accuracy by an FIR filter, and any FIR filter can be approximated by an IIR filter, even though the latter can be much more difficult. The appropriate design criterion obviously depends on the properties of the data and the noise. When it comes to smoothing we usually assume sufficiently oversampled (i.e., smooth) data. If the noise has mainly high frequency components, i.e., if there is little spectral overlap between the data and the noise, we want to maximize stop band attenuation, or minimize stop band energy, while preserving the desired signal as well as possible. In this case we could choose a linear-phase FIR filter designed according to a minimax criterion using the Parks-McClellan algorithm. We could also minimize the stop band energy (i.e., minimize noise power in the stop band) by choosing a least squares method. A mix between the two criteria (minimax and least squares) is possible by choosing a constrained least squares design, which minimizes stop band energy while constraining the maximum approximation error in the pass band.

If the noise spectrum significantly overlaps with the signal spectrum, a more careful approach is required, and brute-force attenuation will not work well because either you leave too much noise (by choosing the cut-off frequency too high) or you distort the desired signal too much. In this case Savitzky-Golay (S-G) filters may be a good choice. The price to be paid is mediocre stopband attenuation, but one advantage is that some signal properties are preserved very well. This has to do with the fact that S-G filters have a flat pass band response, i.e.,

$$\frac{d^kH(e^{j\omega})}{d\omega^k}{\huge |}_{\omega=0}=0\quad k=1,2,\ldots,r\tag{1}$$

where $r$ is the order of the approximating polynomial and $H(e^{j\omega})$ is the filter's frequency response. Property $(1)$ guarantees that the first $r$ moments of the input signal are preserved in the output, which means that the width and height of peaks in the desired signal are well preserved.

Of course there's also a compromise between the two types of smoothing filters discussed above (high stopband attenuation and S-G). We could design an FIR filter with a certain degree of flatness at $\omega=0$ and use the remaining degrees of freedom to maximize stop band attenuation, or minimize stopband energy. In the case of FIR filters, the resulting design problem is sufficiently simple (and convex), and general optimization routines available in several software packages can be used to obtain the optimal filter for the given application.

For the ones interested in the theory of S-G filters, the most relevant references I can recommend are the following:

the original paper by Savitzky and Golay
a tutorial paper by Ronald W. Schafer
section 8.3.5 of Orfanidis' DSP book


How Can I Use MATLAB to Solve a Total Variation Denoising / Deblurring Problem? The Total Variation Denoising Problem is given by:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \lambda \operatorname{TV} \left( x \right) $$

Where $ \operatorname{TV} \left( \cdot \right) $ is the Total Variation Norm.

How could one solve this using MATLAB?

I will solve this for 1D but it could easily generalize into 2D.
The nice thing about the TV Norm that it can be re formulated by the $ {L}_{1} $ norm of the Derivative Operator:
$$ \operatorname{TV} \left( x \right) = \sum_{i = 1}^{N - 1} \left| {x}_{i + 1} - {x}_{i} \right| = {\left\| D x \right\|}_{1} $$
Where $ D $ is the matrix form of the Derivative Operator.
So the whole problem can be formulated as:
$$ \arg \min_{x} \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \lambda \operatorname{TV} \left( x \right) = \arg \min_{x} \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \lambda {\left\| D x \right\|}_{1} $$
For the above there are many solvers based on ADMM or Proximal Gradient Descent.
You could even use the Sub Gradient Method to solve the above as the Sub Gradient is given by:
$$ {A}^{T} \left( A x - y \right) + {D}^{T} \operatorname{sign} \left( D x \right) $$


Recommendation for studying filters I have a project to do about what happens to a periodic function when we pass it through a low-pass, high-pass and band-pass filter.
I have no expressions for the filters or the function, I just have to analyse graphics.
I already concluded that the low-pass filter passes the zones of the graphic that have a low variation, and that the high pass filter passes the zones of the graphic that have a high variation. In the case of the band-pass filter the resulting wave consists on various sinusoids (it just passes certain intervals of variations).
Ok, I now have the general idea, but I was looking for a text that provided me  more rigorous explanation in terms of the frequencies, and their relation to the variation of a function.
Does anyone here knows a book, or a paper that covers this specific topic of signals filtering?
Thanks!

If you have some numerical  methods background,  Richard Hammings, Digital Filters, is a good place to start.  The book is from 1977 but the math hasn’t changed.  The book has also been released by Dover so it has a very low price, unlike Oppenheim and Schaefer’s book and also unlike Oppenheim and Schaefer, is compact.
The math level isn’t advanced.

Hamming is one of the founders of Information Theory.  He essentially invented parity codes. He was also the first president of SIAM.


Integration with infinite/finite limits as a form of summation in DSP notation for discrete signals While reading some DSP literature, I've noticed integrals of the following form.  An example integral (Equation 7.17) is given on pg. 124 of a book that I've found very challenging to read (Seismic Inverse Q filtering):

$u(t) = h(t)\int\limits_{ - \infty }^\infty  {u(\tau ,t)d} \tau$

This is exactly how the author has written the equation.  Gauging the author's intent, I think that $u(\tau ,t)$ is really another function, so in the above, let's re-write the equation as the following:

$u(t) = h(t)\int\limits_{ - \infty }^\infty  {p(\tau ,t)d} \tau$

In the above, $u(t) \neq p(\tau ,t)$, so I think that this clears up some of the notation.

This appears to be notation associated with a continuous signal $p(\tau, t)$ that is of infinite length.  For a discrete signal $p[\tau, t]$ with a finite length $N$, how do I evaluate the integral?

Case #1

Is the integration treated as a summation?

$u[t] = h[t]\sum\limits_{i = 1}^N {p[i,t]}$

Case #2

Is the integration treated as a "classical" integration using the trapezoidal rule?

How do I know which of these cases apply for the discrete signal case?  Is this really an issue of notation?

Another Example

Another example can be found on pg. 128 of the same book (Equation 7.21):

$\tilde U(\tau ,\omega ) = U(\tau ,\omega )\frac{1}{{\Lambda (\tau ,\omega )}}\exp \left[ {i\int\limits_0^\tau  {\left( {{{\left( {\frac{\omega }{{{\omega _h}}}} \right)}^{\gamma (\tau ')}} - 1} \right)} \omega d\tau '} \right]$

When the integration is being done over $\tau'$, is this integration a summation or a "classical" integration?

Yet another example

On pg. 9 of  "R. Wang, Introduction to Orthogonal Transforms: With Applications in Data Processing and Analysis. Cambridge University Press, 2012", the author states that the energy in a continuous signal $x(t)$ is given by:

$E = \int\limits_{ - \infty }^\infty  {{{\left| {x(t)} \right|}^2}} dt$

The energy in a discrete signal $x[n]$ is:

$E = \sum\limits_{n =  - \infty }^\infty  {{{\left| {x[n]} \right|}^2}}$

But shouldn't a discrete signal have finite lower and upper bounds of summation?  What is really meant here using infinite lower and upper bounds?  Shouldn't the trapezoidal rule be used here?

So here is the gist of my question: Using this notation, is integration done using the trapezoidal rule for discrete signals, or is it simply a form of summation?  What is the general use of this type of notation in DSP?

Given an equation for a continuous function integration, how does this relate to integration in the context of a discrete function?

In the absence of any further information, we must assume that tradition plays a role in the numerical implementation.  The language of mathematics implies that an integration is an integration, so we can assume that the trapezoidal rule holds:

$\int_a^b f (x){\mkern 1mu} dx \approx \frac{1}{2}\sum\limits_{k = 1}^{N} {\left( {{x_{k + 1}} - {x_k}} \right)} \left( {f({x_{k + 1}}) + f({x_k})} \right)$

However, the intent of the author is another mitigating factor, so careful reading is required.  There are often some pathological issues of notation, but these are strange.

After some numerical experimentation, I've found that in general, an integration is an integration, no matter which field of inquiry.  The language of mathematics is universal.


ARMA vs. AR and then what? Sorry if this sounds elementary but I am struggling to grasp the physical idea behind ARMA (auto-regressive, moving average) process. The "AR" part is intuitive and so is "MA", but put together?

If I model my time series using "AR", I can predict the next sample using linear terms of the previous samples and the difference between the predicted and the real one should then have a normal distribution provided the noise was guassian white (the unpredictable).
But what does it mean to say that my time series can be ARMA modelled?

thanks

The arma approach is to model the current output of the system as the sum of past outputs and past inputs explicitly. The assumption of gaussian model for the noise statistics still can be used for the unpredictable signal which cannot be modeled as arma.

From system's frequency response spectrum modeling aspect: an ar model is able to model only the spectrum peaks as it fits the peaks better (using poles) and the ma model is able to model the valleys better (using zeros) in the spectrum. Hence, by using an arma model we will be better enable both peaks and valleys using a reasonable order for both ar and ma.

However, it is always possible to model a zero (in z-domain) using a infinite number of poles so an approximation to an arma model (whose parameter estimation is relatively involved) is to use a high order ar model (whose paramter estimation is relatively easy).


Plotting Frequency Response Magnitude and Phase for first order all pass filter Having trouble plotting the frequency response characteristics for first order all pass filter. The Magnitude is expected to be constant across entire freq and Phase is expected to be only decreasing (with jump from -pi to pi). I don't see either of this happening, provided plot from octave.
$H_{ap}(z) =  {z^{-1} - a^* \over 1 - az^{-1}}\tag{1}$
Above is condition to enable all pass filter.
code.m
num = [0.8,1]
den = [1;0.8]
TF = tf(num,den)
[zz,pp,kk] = tf2zp(TF)
figure(1);
bode(TF);
figure(2);
pzmap(TF);

Octave output:
Octave output:

Transfer function 'TF' from input 'u1' to output ...

      0.8 s + 1
 y1:  ---------
       s + 0.8

Continuous-time model.
zz = -1.2500
pp = -0.80000
kk =  0.80000

As pointed out by Hilmar in a comment, you should remain in the discrete domain. This is most easily done by using the functions freqz() and zplane():

b = [.8,1]; a = [1,.8];
[H,w] = freqz(b,a,2048);
subplot(2,1,1), plot(w/pi,abs(H));
subplot(2,1,2), plot(w/pi,angle(H));
figure, zplane(b,a)

Note that for a first-order all-pass filter it's also quite straightforward to come up with analytical solutions for the magnitude and phase of the frequency response.


Count Matches to a Kernel So, I have this problem where I want to apply a kernel to an image and count the number of matches that happened.

So for example, if I have the kernel:
$$\begin{bmatrix}
1 & 2 & 1\\
1 & 2 & 1\\
1 & 2 & 1\\
\end{bmatrix}$$
And the image:
$$\begin{bmatrix}
1 & 1 & 2 & 1 & 2\\
1 & 1 & 2 & 1 & 2\\
1 & 1 & 2 & 1 & 2\\
\end{bmatrix}$$
The result should be (ignoring the borders):
$$\begin{bmatrix}
3 & 9 & 0\\
\end{bmatrix}$$
Explanation: when the kernel is in the first position (2,2), the 3 left ones match to the kernel, in the middle position (3,2) all 9 elements match and in the right (4,2) there is no match.

I ended up doing two convolutions (or correlations), one to count the numner of ones and one to count the number of twos.

For example:
First convolution (to count where value is 1):
$$\begin{bmatrix}
1 & 0 & 1\\
1 & 0 & 1\\
1 & 0 & 1\\
\end{bmatrix}$$
$$\begin{bmatrix}
1 & 1 & 0 & 1 & 0\\
1 & 1 & 0 & 1 & 0\\
1 & 1 & 0 & 1 & 0\\
\end{bmatrix}$$
Result:
$$\begin{bmatrix}
3 & 6 & 0
\end{bmatrix}$$
Second convolution (to count where value is 2):
$$\begin{bmatrix}
0 & 1 & 0\\
0 & 1 & 0\\
0 & 1 & 0\\
\end{bmatrix}$$
$$\begin{bmatrix}
0 & 0 & 1 & 0 & 1\\
0 & 0 & 1 & 0 & 1\\
0 & 0 & 1 & 0 & 1\\
\end{bmatrix}$$
Result:
$$\begin{bmatrix}
0 & 3 & 0
\end{bmatrix}$$

At the end I add up the results to get the expected result. The problem is that if my image has 256 values, I have to do that 256 times and it becomes 256x slower than a convolution. Still, this ends up being faster than doing for loops and checking all kernel positions. This is partly due to the fourier transform that makes convolution faster.

My question is if it is possible to apply this operation faster? Maybe by just modifying the convolution a little bit?

EDIT:
My kernel has maximum size 64x64 and my image has size around 300x300.
If I use the standard kernel convolution (Direct Convolution), then it will try to match the kernel for every position of the image. To match the kernel it takes 64x64 comparisons, and it will do that for every position in the image (300x300), so this will take 300x300x64x64 comparisons (total 368.640.000 operations).

If I compute first the fourier transform of the images (FFT Convolution), then I can do the multiplication directly and it will cost only (300x300 operations). There is a cost to compute the fourier of the images in the first place, so at the end it will take around 300x300xlog(300)xlog(300) operations (this takes a total of around 500.000 operations).

OpenCV has a convolution method (filter2D) and a template matching method (matchTemplate) that are very fast, and I believe they use fourier plus some other optimizations.

The function matchTemplate in OpenCv has different methods to compute the match https://www.docs.opencv.org/2.4/doc/tutorials/imgproc/histograms/template_matching/template_matching.html, but none of them work for what I need.

What I need is a function to count the matches, and I want to know if it is possible to that similar to other template matching methods in opencv (that I assume use the FFT Convolution).

An approximate result should be fine, so if the output for the first example is, lets say, [3.28, 9.45, 0.32], it shouldn't be a problem.

Besides what is discussed below (see for instance the above Cris Luengo's comment): for reasonably-sized images, and modern computing techniques, I can see no reason in practice why the basic template match counting could be slow, and why one should look for faster, imprecise and incomplete alternatives.

[UPDATE 20180730]
You are trying to compare matches in a pair of matrices. Ontologically, the comparison operation  is not a linear operation. Hence, no convolution operation, naturally linear, seems be the easiest approach. Neither the FFT or Fourier stuff in general, as they are quite related to linearity.

IMHO, you ought to look at the concept of Template Matching, and avoid the notion of kernel, which is a debated notion.

Instead, a simple counting approach could do the job. For instance in Matlab:

%% The OP kernel and image
% dataKernel = [1 2 1;1 2 1;1 2 1;];
% dataImage = [1 1 2 1 2;1 1 2 1 2;1 1 2 1 2;];

%% A random kernel and image
dataKernel = floor(256*rand(64,64));
dataImage = floor(256*rand(300,300));

%% Compute sizes
[nRowKernel,nColKernel] = size(dataKernel);
[nRowImage,nColImage] = size(dataImage);

[nRowCrop] = nRowImage - nRowKernel+1;
[nColCrop] = nColImage - nColKernel+1;

lMatch = zeros(nRowCrop,nColCrop);
dataImageCrop = zeros(size(dataKernel));

tic;
for iRowCrop = 1:nRowCrop;
    for iColCrop = 1:nColCrop;
        dataImageCrop = dataImage(iRowCrop:iRowCrop+nRowKernel-1,iColCrop:iColCrop+nColKernel-1); % Crop a portion of the image
        lMatch(iRowCrop,iColCrop) = length(find(dataImageCrop - dataKernel == 0));
    end
end

Here, I just used an ad-hoc displacement for the crop, but than could be automated easily. No product, just differences. If your data are limited size integer, there might be more options for speed. But even in Matlab, this is relatively fast.

FFT-based cosine similarity for fast template matching seems to better explain what is alluded too in OpenCV performance on template matching (provided by the OP). If the image is of size $I_x \times I_y$, and the template is  $T_x \times T_y$, with $T_x < I_x$ and  $T_y < I_y$, you can expect a number of operations of:

$O(I_xI_y  T_x T_y)$ with the standard comparison,
$O(I_xI_y\log I_x \log I_y)$ with a Fourier method.

From afar, Fourier might seem to be beneficial whenever $\log I_x \log I_y <  T_x T_y$. But:

this is only an asymptotic result the constant behind the $O(\cdot)$ can make a huge difference
then you have to find maxima, which add to the overall cost
the hardware implementation (pipeline, CPU/GPU) can made a lot of difference in practice, compared to number of operations
indeed, the first one only require adds/subtracts and absolute values on integers,  the second complex operations on floats; thus, there might be huge variations between implementations
a difference of pixels, like the sum of absolute difference (SAD) does not translate easily in Fourier because it is orthogonal
implementing fast correlations, the FFT might be able to provide you with a "best match", but it is not local, so I cannot imagine right now how this can directly give you a count by patch location.

So far, to me:

if you only want a best match, then you can do Fourier on the whole image;
if you really want (even approximate) counts at all template locations, I do not see an easy and useful Fourier implementation for that so far.

However, you can speed the matching up, with several tricks, like using the symmetry of the template (if any), the integral image, multiscale pyramid-like matching, etc. A few additional pointers:

Fast Template Matching
Fast Template Matching With Polynomials
FasT-Match: Fast Affine Template Matching - The Computer Vision ...
Very Fast Template Matching


Text in image is low or high frequency component? Text present in captured images is high frequency or low frequency component?Why?

The edges of the words will probably contribute to a high frequency, the background of the text and the portion within the boundaries of the letters of the text will be otherwise.


RMS value of a signal with a DC offset in MATLAB t = 0:0.001:1;
y1 = 400 * sin(2*pi*20*t);
y2 = 200 + 400 * sin(2*pi*20*t);
plot(t,y1)
hold on
plot(t,y2,'-g')
grid on;
rms_y1 = sqrt(mean(y1.^2))
rms_y2 = sqrt(mean(y2.^2))

In the above MATLAB code there are two signals y1 and y2.

y1 is a sinusoidal signal 20 Hz with 400V amplitude.

y2 is a signal which is nothing but y1 + 200V DC offset.

So I was thinking the "rms of y2" would be equal to "rms of y1 + 200". Since rms of a constant is the constant itself.

But rms of y1 is 282.7014; and rms of y2 turns out to be 346.2948 (not 282 + 200).

What am I doing wrong here?

RMS values don't simply add up. The RMS value of $N$ discrete samples $x[n]$ (assuming $0\le n<N$) is

$$x_{RMS}=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}x^2[n]}\tag{1}$$

Now if you have $y[n]=x[n]+c$ with some constant $c$, the RMS value of $y[n]$ is

$$\begin{align}y_{RMS}&=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}y^2[n]}\\&=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}(x[n]+c)^2}\\&=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}(x^2[n]+2cx[n]+c^2)}\\&=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}x^2[n]+\frac{2c}{N}\sum_{n=0}^{N-1}x[n]+\frac{1}{N}\sum_{n=0}^{N-1}c^2}\tag{2}\end{align}$$

If we define the mean value of $x[n]$ as

$$m_x=\frac{1}{N}\sum_{n=0}^{N-1}x[n]\tag{3}$$

then from $(2)$ the square of the RMS value of $y[n]$ can be expressed as

$$y^2_{RMS}=x^2_{RMS}+2cm_x+c^2\tag{4}$$

So if you assume that the mean of $x[n]$ equals zero, the squares of the RMS values add up:

$$y^2_{RMS}=x^2_{RMS}+c^2\tag{5}$$

You can easily verify Eq. $(5)$ using your example:

$$(282.7014)^2+(200)^2=(346.2948)^2\tag{6}$$


Unable to understand the derivation of the update equation for LMS I am trying to follow the derivation of the Least Mean Square https://en.wikipedia.org/wiki/Least_mean_squares_filter#Proof
 but I cannot get the update rule. I am stuck in the following steps and shall appreciate help where I am going wrong.
The transmitter generates a sequence of random input data $\{s_n\}$, each element of which belongs to an alphabet $\mathcal{A}$. The data sequence is sent through the channel whose output $y_n$ is observed by the receiver. The input/ output relationship of the channel can be written as:
\begin{align}
y_n &= \sum_{k=1}^L h_k s_{n-k} + v_n,\nonumber\\
&= \mathbf{h}^H \mathbf{s}_n + v_n. \tag{1}
\end{align}
where the additive noise $v_n$ is assumed to be stationary, Gaussian, and independent of the
channel input $s_n$. We also assume that the channel is stationary, moving-average and has finite length. The output of the equalizer is:
\begin{align}
z_n &= \sum_{k=1}^L {w}_k y_{n-k},\nonumber\\
&=\mathbf{{w}}^H \mathbf{y}_n = \hat{s}_n \tag{2}
\end{align}
where the superscript $H$ denotes conjugate transpose, $\mathbf{y}_n =[y_n,y_{n-1},\ldots,y_{n-L+1}]^T$ is the signal to be equalized, $z_n$ is the output of the equalizer and $\mathbf{{w}} = [{w}_1,{w}_2,\ldots,{w}_L]^T$ is the vector of equalizer coefficients.

The cost function for LMS is $J = E[|e_n|^2]$ which is minimized.
\begin{align}
e_n &= s_n - \hat{s}_n,\\
&=s_n - \mathbf{{w}}^H \mathbf{y}_n \tag{3}
\end{align}
We need to take the gradient of the cost function and then equate it to zero to yield the equation for the equalizer coefficients. The gradient is
\begin{align}
\nabla J(\mathbf{w}) &= \nabla E[e_ne^*_n],\\
&=2E[\nabla(e_n)e^*_n]
\tag{4}
\end{align}
Now, $\nabla(e_n) = \nabla(s_n - \mathbf{{w}}^H \mathbf{y}_n) = -\mathbf{y}_n$. Substituting this into Eq(5) and equating to zero we get,
\begin{align}
\nabla J(\mathbf{w}) &= \nabla E[e_ne^*_n],\\
&=2E[-\mathbf{y}_ne^*_n]
\tag{5}
\end{align}
Using Eq(3),
\begin{align}
\nabla J(\mathbf{w}) =0,\\
2E[-\mathbf{y}_ne^*_n] =0,\\
\mathbf{y}_ne^*_n = 0,\\
\mathbf{y}_n(s_n^* - \mathbf{{w}} \mathbf{y}_n^*) =0,\\
\mathbf{{w}} = {(\mathbf{y}_n \mathbf{y}_n^*)}^{-1}(\mathbf{y}_n s_n^*)
\tag{6}
\end{align}

I cannot understand how to plug in the $\mu$, the step size. This Eq(5) is completely different from the update equation :$\mathbf{{w}}_{n+1}=\mathbf{{w}}_n + \mu e_ns_{n}^*$

so with an LMS filter, we have a time-variant $N$-tap FIR filter:

$$ y[n] = \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] $$

$x[n]$ is the input signal, $y[n]$ is the FIR output, and $h_n[k]$ are the FIR tap coefficients at the time of sample $n$.

with an LMS filter, we also have another input called the desired signal: $d[n]$.  we want our LMS filter to adapt so that the output $y[n]$ follows the desired signal $d[n]$.  from those two signals, we define an LMS error signal:

$$ e[n] \triangleq y[n] - d[n] $$

and this Least Mean Square adaptive filter will adjust its tap coefficients, $h_n[k]$, in such a manner as to eventually minimize the square of this error signal $e[n]$.

so the error signal is:

$$\begin{align}
 e[n] &= y[n] - d[n] \\
      &= \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \\
\end{align}$$

and the square of this error signal (which we want to lessen in an LMS filter) is:

$$ \big(e[n]\big)^2 = \left( \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \right)^2 $$

now, one thing Bernard Widrow must have been thinking is that, since the inputs $x[n]$ and $d[n]$ are given, the only controls we can adjust to lessen this square error $e^2[n]$ are the $N$ tap coefficients $h_n[k]$ for $0 \le k \le N-1$.  with $x[n]$ and $d[n]$ fixed, then $e^2[n]$ is a function, $f(\cdot)$, of the $N$ tap coefficients $h_n[k]$.  so at time $n$ this function of $N$ variables looks like

$$\begin{align}
 e^2[n] & \triangleq f\big(h_n[0], h_n[1], h_n[2], ... h_n[N-1] \big) \\
        &= \left( \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \right)^2 \\
\end{align}$$

at time $n$, the derivative of $f(\cdot)$ with respect to the $m$-th tap coefficient is

$$\begin{align}
 \frac{\partial}{\partial h_n[m]} \, e^2[n] &= \frac{\partial}{\partial h_n[m]} \, f\big(h_n[0], h_n[1], h_n[2], ... h_n[N-1] \big) \\
    &= \frac{\partial}{\partial h_n[m]} \, \left( \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \right)^2 \\
    &= 2 \, \left( \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \right) \cdot x[n-m] \\
    &= 2 \, e[n] \cdot x[n-m] \\
\end{align}$$

now from this we can set all of these derivatives to zero and set up $N$ equations and $N$ unknowns and try to solve the mess.  i am sure Widrow did this.  but it doesn't lead to a very cheap (in terms of computational cost) algorithm.

but for each $m$, where $0 \le m \le N-1$, the quantity

$$ \frac{\partial}{\partial h_n[m]} \, f\big(h_n[0], h_n[1], h_n[2], ... h_n[N-1] \big) = 2 e[n] x[n-m] $$

is one component of the vector which is the gradient of $f(\cdot)$.

for the moment, try to visualize this as a function of two variables: $f\big(h_n[0], h_n[1]\big)$.  this is a surface above the point in a plane with coordinates $(h_n[0], h_n[1])$.  the gradient

$$\begin{align}
 \left(\tfrac{\partial f(\cdot)}{\partial h_n[0]}, \ \tfrac{\partial f(\cdot)}{\partial h_n[1]}\right) &= \Big(2 e[n] x[n-0], \ 2 e[n] x[n-1] \Big) \\
&= 2 \, e[n] \Big(x[n], \ x[n-1] \Big) \\
\end{align} $$

is the 2-dim vector that points in the direction of the steepest increasing slope of this surface.  so if $\tfrac{\mu}{2}$ is a tiny positive number, moving a tiny bit in the direction of the gradient should increase this function a tiny bit.

$$ f\big(h_n[0], \ h_n[1]\big) < f\big(h_n[0] + 2 e[n] x[n] \tfrac{\mu}{2}, \ h_n[1] + 2 e[n] x[n-1] \tfrac{\mu}{2} \big) $$

but $f(\cdot)$ is the square of error and we want to decrease it, so we move in the opposite direction of the gradient:

$$ f\big(h_n[0], \ h_n[1]\big) > f\big(h_n[0] - 2 e[n] x[n] \tfrac{\mu}{2}, \ h_n[1] - 2 e[n] x[n-1] \tfrac{\mu}{2} \big) $$

that's with two taps ($N=2$).  extending the idea of this "surface" to more dimensions than two, we expect

$$ f\big(h_n[0]-e[n]x[n]\mu, ... ,h_n[k]-e[n]x[n-k]\mu, ... \big) < f\big(h_n[0], ... ,h_n[k], ... \big) $$

so, for a small stepsize $\mu$, the adaptation algorithm (or coefficient "update equation") is

$$ h_{n+1}[k] \ \leftarrow \ h_n[k] \ - \ (\mu \, e[n]) \, x[n-k] $$

that is the simple unnormalized LMS adaptive filter.  you want $\mu$, also known as the adaptation gain, to be large enough so that the LMS filter adapts at an acceptably quick rate, but small enough so that it doesn't overshoot in the adaptation and become unstable.

because the stepsize is really $(\mu \, e[n]) x[n-k]$ and both $e[n]$ and $x[n-k]$ are expected to be proportional in magnitude to the mean $\big|x[n]\big|$, this means louder signals will cause the adaptation to happen faster than quieter signals.  so from this, the Normalized LMS adaptive filter is defined to divide that term by some measure of the mean of $\big|x[n]\big|^2$.

so the Normalized LMS adaptation algorithm (or coefficient update) is

$$ h_{n+1}[k] \ \leftarrow \ h_n[k] \ - \ \frac{\mu \, e[n]}{\overline{x^2[n]} + \epsilon} \, x[n-k] $$

where $\overline{x^2[n]}$ is some measure of mean of $x^2[n]$ for the samples most recent to time $n$.  you can simply square $x[n]$ and run that through a simple first-order IIR filter to get a measure of the mean:

$$ \overline{x^2[n]} \ \leftarrow \ p \, \overline{x^2[n-1]} \ + \ (1-p) x^2[n] $$

where $p$ is the pole of the one-pole IIR filter and $0 < (1-p) \ll 1$.  the quantity $\epsilon > 0$ is a tiny amount to add to the denominator so that we never divide by zero, even if the input $x[n]$ goes to zero.


Comparison between Fourier transform, short-time Fourier transform and wavelets What is the difference between the Fourier transform, short-time Fourier transform and wavelets?

All three transforms are inner product transforms, meaning the output is the inner product of a family of basis functions with a signal. The parametrization and form of the basis functions determine the properties of the transforms.The number of basis functions for a complete picture (i.e. a result that contains enough information to reconstruct the original signal) is  the same for all three cases. That is a result of the Heisenberg uncertainty or a bit more general, conservation of information.

For the standard Fourier transform the basis functions are simply the complex oscillations
$$b_\omega(t) := \exp(i \omega t)$$
where $t$ is the time axis of the signal and $\omega$ is the single frequency parameter that determines the basis function in the family. There is one basis function for every $\omega$. The fourier transform of the signal $s(t)$ is then simply the inner product written as an integral:
$$\mathcal{F}\left\{s(t)\right\}(\omega)=\langle b(\omega,t),s(t)\rangle=\int_{-\infty}^\infty \exp(-i \omega \tau) s(\tau)d\tau$$
The negative sign in the exponential comes from the complex conjugation in the general inner product definition $$\langle a(t),b(t)\rangle=\int_{-\infty}^\infty a(\tau)^*b(\tau) d\tau$$

The short time fourier transform adds a time dimension to the base function parameters by multiplying the infinitely long complex exponential with a window to localize it. The base functions are then something like $$b_{(\omega,t_0)}(t):= w(t-t_0)\exp(i \omega t)$$ where $w(t)$ is the window functions that vanishes outside some interval and $(\omega,t_0)$ are the time-frequency coordinates of the base function in the family. The inner product becomes then
$$\mathcal{S}\left\{s(t)\right\}(\omega,t_0) = \int_{-\infty}^\infty w(\tau)^* \exp(-i\omega \tau) s(\tau)d\tau$$
and gives you time and frequency information. The actual tradeoff between time and frequency is determined by the choice of the window function. The result of this transform can also be regarded as a filterbank with bandpass filters that have the fourier transform as the window $w(t)$ as frequency response, but shifted to the center frequency $\omega$. All the filters therefore have the same bandwidth.

The wavelet transform does away with the constant bandwidth constraint and adapts the window size to the frequency. And that happens in a specific scale invariant way that doesn't even need the complex modulation anymore. You start with a generic base function that is localized and oscillates, and it also has zero mean (i.e. the integral over the complete space is 0). That is your wavelet. Now if you scale the wavelet in time direction the oscillation frequency changes as well, so you can control localization and oscillation with a single parameter that links the two. The family of base functions is therefore something like $$b_{(\sigma,t_0)}(t) = w\left(\frac{t-t_0}{\sigma}\right)$$ where $w$ is the original wavelet and $\sigma$ the scale parameter. The inner product becomes then $$\mathcal{W}\left\{s(t)\right\}(\sigma,t_0)=\int_{-\infty}^\infty w\left(\frac{t-t_0}{\sigma}\right)^* s(t)$$
Again this gives you a 2-dimensionally parametrized result that can also be seen as a filterbank, but this time the bandwidth of the filters is proportional to the center frequency of the bands. The time-frequency plane is therefore partitioned non-uniformly and more suitable to many real-world signals where the typical time constants of processes is naturally scale dependent.


What does the frequency band mean when it comes to finding aliases? The time signal which i'm trying to find the aliases for is:

$$x:{\mathbb R}\rightarrow {\mathbb R}\\\   x(t)=\cos(50t) +2\cos(70t).$$

If the sample period is $T_s = \frac{\pi}{60}$ then according to Nyquist -Shannon sampling theorem (which btw my professor failed to prove) there is/are a signal(s) which after sampling will be equal to the sampled version of the above signal, if we sample those with the same sample frequency in the frequency band $[-55, 55]$.

I don't understand the meaning of the last sentence , what does a frequency band means here?

Let's use a simple example: A sine wave of 400 Hz. The Fourier spectrum has two components: one at 400 Hz and one at -400 Hz (since it's a real valued time signal).

Now let's sample this at 1 kHz. Sampling in one domain is periodic repetition in the other domain. So the sampled signal consists of your two original frequencies, -400 Hz and +400 Hz repeated again and again in 1 kHz intervals. So you get components at +-400 Hz, +-600 Hz, +-1400 Hz, +-1600 Hz, etc.

Any of these mirror frequencies (600, 1400, 1600, 2400, 2600, ...) will give you the exact same samples as the original 400 Hz.


What are some alternatives to Spectrogram charts? Spectrogram renders of audio are limited in their representation of audio. For example, one could not use a Spectrogram to re-render the original audio.

What are some alternatives to Spectrogram charts for rendering representations of audio?

Preferably this would be methods which are more accurate (in every possible way/in every possible audio dimension) for representing audio in graphical format (an image)?

A good measure of accuracy would be for example;

Could the image be used to (at least to some extent) re-create the audio file?
Could the image be used to match the similarity between two songs in some or all dimensions?

And, are there any scripts/tools/software's that can be used for converting audio into an image format?

The ltfatpy 1.0.16 package

  is a partial Python port of the Large Time/Frequency Analysis Toolbox
  (LTFAT), a MATLAB®/Octave toolbox for working with time-frequency
  analysis and synthesis.

Among linear and quadratic time-frequency methods, there is a large number of options for sharper analysis tools converting a 1D signal into 2D data. You can even get more information by decomposing onto 3D or 4D spaces (coming back on that later). Quadratic (or higher-order) representations are usually more difficult to invert than linear ones.

And here is a documentation for ltfatpy.


Noise Equivalent Bandwith of IIR and FIR Filters? How do you calculate the noise equivalent bandwidth for an IIR or FIR filter?

Hm, isn't it the same as for any filter?

$$
B_\text{eq} = \frac{1}{|H(e^{j\omega_0})|^2} \int_0^\infty |H(e^{j\omega})|^2 d\omega
$$

where $H(e^{j\omega})$ is the frequency response, and $\omega_0$ is the max-abs frequency. For FIR filters, you can just use e.g. freqz in MATLAB/Octave. For IIR filters, you'd have to do analysis or take measurements (or see if you can get away with truncating)


Level Correction For Digital Down-Convertion? I am trying to understand a subtle (sorry if this is obvious) difference in the computed output power after a DDC process. Here is the story:

Suppose I generated a 10 MHz signal from a microwave source at -10 dBm and feed it to the a digitizer. This corresponds to a peak-to-peak voltage of 200 mV as measured by an oscilloscope for 50 Ohm lines and termination.

With the digitizer and by performing DDC, I am getting around +4 dBm as a result, 14 dBm above the expected -10 dBm of the source. Why is that?

I tried to go through the math step-by-step:

With the 200 mV input peak-to-peak signal, by multiplying this signal with a cos (for I) and -sin (for Q) with amplitude equals 1, I will obtain (suppose the signal has zero phase) I_DC = 50 mV and Q_DC = 0 mV after filtering out the high frequency part, just by trigonometry.

The power is then computed as P=I^2+Q^2= 2.5 mW. Taking log, we get P_log = 10*log(2.5mW) ~ -26 dB = 4 dBm. This is around 14 dBm higher than what I originally generated (-10 dBm).

If I wasn't missing something in the calculation, it seems I need a level correct factor of 14 dBm. But where did this gain come from?
I guess by duplicating the signal to generate the I and Q there is a 6 dB gain. Then where is the remaining 8 dB from?

Appreciated if you can shine some light on it.

Edit1: (Answer?)
Following the first reply, I might need to consider the 50 Ohm, so it will bring the power down by a factor 10*log(1/50), and it gives now P_dBm = 10*log(2.5e-3*1000/50) = -13.0103 dBm. I am now 3.0103 dB off, and it seems it comes from the trigonometry factor of half. By putting it back, I come back to the -10 dBm as desired.

I think you miscalculated the power generated by your signal source; note that you need to take the load impedance into account:

$$
\begin{align}
P = \frac{V_P^2}{2R} &= \frac{(100\ \text{mV})^2}{100\ \Omega} \\
&= \frac{0.01\ \text{V}^2}{100\ \Omega} \\
&= 0.1\ \text{mW} \\
\\
\end{align}
$$

$$
\begin{align}
P_{\text{dBm}} &= 10 \log_{10}\left(\frac{P}{1\ \text{mW}}\right) \\
&= 10 \log_{10}\left(\frac{0.1\ \text{mW}}{1\ \text{mW}}\right) \\
&= -10.0\text{ dBm}
\end{align}
$$

Likewise, if your DDC output signal has a peak amplitude of $50\text{ mV}$, then you can repeat the above calculations to yield a power level of $-13\text{ dBm}$.


Format location from LPC magnitude I am confusing when I find the lpc roots magnitude, f1=0.9313, f2=0.9308, f3=0.91 and f4=0.95 but when I look at LPC envelope the f1 is higher than f4 in term of magnitude in dB.
Could anybody explain it ?

I have wrote following MATLAB code :

frame=wave((i-1)*frame_length+1:(i*frame_length));
a=lpc(frame,30);
r=roots(a);
r=r(imag(r)>0.01);           % only look for roots >0Hz up to fs/2
for k=1:length(r)
z(k) = abs(r(k));
ffreq(k)=(atan2(imag(r(k)),real(r(k)))*fs/(2*pi));% convert to Hz and sort
end
[ff,index]=sort(ffreq);
fprintf('\n Start of frame %d \n ',i);
for j=1:length(ff)
fprintf('Formant %d Frequency %.1f --- %.2f\n',j,ff(j),z(index(j)));
end
fprintf('\n End of frame %d \n ',i);

The magnitude of the roots does not give you the format gain. You need to evaluate the frequency response of the filter $H(z)=\frac{1}{A(z)}$, where $A(z)$ is the polynomial returned by the LPC function, at the point on the unit circle corresponding to the desired format frequency. Like this freqz(1,a, angle(r)) for instance.


Apply a Gabor filter to an input image I tried to apply a Gabor filter with a specific scale (according to my values of lambda and sigma, so it is (7x7) and for 4 orientations (0, $\frac{\pi}{4}$, $\frac{\pi}{2}$ and $\frac{3\pi}{4}$) to a input gray scale image.

In my code, three steps are achieved:

Create a Gabor Filter
Read an RGB image, then convert it to gray scale and finaly to double.
Apply the created gabor to the input image (here, i am not sure if my code is true, that's why i need your opinion )

1) -------------- create the Gabor filter (size=7x7 and 4 orientations)

%define the five parameters
theta=....; %either 0 or pi/4 or pi/2 or 3pi/4
lambda=3.5;
gamma=0.3;
sigma=2.8;
psi=0;

sigma_x = sigma;
sigma_y = sigma/gamma;

nstds = 5;
xmax = max(abs(nstds*sigma_x*cos(theta)),abs(nstds*sigma_y*sin(theta)));
xmax = ceil(max(1,xmax));
ymax = max(abs(nstds*sigma_x*sin(theta)),abs(nstds*sigma_y*cos(theta)));
ymax = ceil(max(1,ymax));
xmin = -xmax; ymin = -ymax;
[x,y] = meshgrid(xmin:xmax,ymin:ymax);

x_theta=x*cos(theta)+y*sin(theta);
y_theta=-x*sin(theta)+y*cos(theta);

gb= exp(-.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+psi);

figure(2);
imshow(gb);
title('theta=...');
%imagesc(gb);
%colormap(gray);
%title('theta=...');

2)------------Read the input image

I=imread('piano.jpg');
image_resize=imresize(I, [160,160]);
image_gray=rgb2gray(image_resize);
image_double=im2double(image_gray);
figure(1);
imshow(image_double);

3) -----apply the created above gabor to the input image (recall that i am not sure if the code in this step is 100% true, that's why i need your opinion and your help if you have the correct answer.)

figure(3);
filtered = conv2(image_double,gb);
imagesc(filtered);
colormap(gray);
title('theta=....');

Your code is correct and the results are consistent. You may be surprised by them due to some 'hidden features'.

First,  conv2 returns by default the full convolution, such that the result is the size of the image plus a border of half the size of the kernel (that is, the total size if the size of the image plus that of the kernel). When you interpret your results, be aware of it!

Second, the results represent coefficients that are stronger for a higher correlation between your kernel and your local image patch: as expected you also extract the borders of the image. See in particular your leftmost result showing strong vertical line.

Last, imagesc scales by default the scale between the highest to the lowest coefficient. That's why in the leftmost result, you mainly see the border.

There are different options to conv2described in help conv2 which allow to control this behavior.

Be aware also that there many different definitions of kernels for detecting edges, such as log-Gabors

.

If you are interested in a full implementation (in python) you may have a look at:  https://pythonhosted.org/LogGabor/ (shameless self-plug ).


Improve spectrogram appearance using windowing functions MATLAB I have a program which plots data obtained in a NI USB DAQ realtime to a spectrogram plot using the surf function. I'd like to improve the appearance of the plot using some sort of windowing function. My questions are: Which windowing function is the best to use for this? and How do I implement it?
Here's a sample of my current code:

fprintf('Initializing DAQ\n');
s = daq.createSession ('ni');
addAnalogInputChannel(s, 'Dev1','ai1', 'Voltage');

s.NumberOfScans = Fs*Record_Time;
s.Rate = Fs;
s.Channels.Range = [-1 1];
s.DurationInSeconds = 5;

Spect_nFFT = 2^8;
num_overlap = 19*round(Spect_nFFT/20);
start_sample = 1;

DataRdy = true;
% Run in background session

Lh = s.addlistener('DataAvailable', @(src,event) CaptureData(event.TimeStamps, event.Data));
s.NotifyWhenDataAvailableExceeds = 100000*.1;
s.IsContinuous = true;
s.startBackground();
NFFT = 2^nextpow2(L);
P = zeros(8193,50);
F = Fs/2*linspace(0,1,NFFT/2+1);
g=0;

while g==0

    while DataRdy == 0
        pause(.01)
    end
    pause(.5)

    Min_Threshold = 10^12;

    start_sample = 1;
    NFFT = 2^nextpow2(L); % Next power of 2 from length of myRecording
    FFT_Out = fft(myRecording(start_sample:end),NFFT)/size(myRecording(start_sample:end),1);
    FFT_Out = FFT_Out(1:NFFT/2+1);

    A = sqrt(FFT_Out.*conj(FFT_Out));
    P = [P(:,2:end), sqrt(FFT_Out.*conj(FFT_Out))];
    P(P<=max(max(P))/Min_Threshold)=max(max(P))/Min_Threshold;
    surf(1:size(P,2),F,10*log10(P),'edgecolor','none'); axis tight;

    % Create updating spectrogram
    view(0,90);

    xlabel('Time');
    ylabel('Hz');

    drawnow;

    DataRdy = 0;
    uicontrol('Style', 'pushbutton', 'String', 'Exit',...
    'Position', [20 20 50 20],...
    'Callback', 'g=g+1');

end
delete(Lh);
stop(s);

Thanks for any help!

Windowing is nothing more than element-wise multiplication of your signal by window function.

Let's assume that you want to apply Hanning window. In your case signal is stored inside myRecording(start_sample:end).

Create the window values vector of your signal length: win =
hann(length(myRecording(start_sample:end)))
Multiply by your signal and store it in variable: sig = win.*size(myRecording(start_sample:end).
Perform the FFT, but keep in mind that now you cannot simply normalize but your signal length (as you do within FFT function), but you should (very simply speaking) divide by sum of your window values: sum(win). Probably you will notice that in case of no window applied = rectangular window, that is also true.

Answering second part of your question, which window to choose? You can start from this great white paper (on figure 8 you have recomendations for different window types): Windowing: Optimizing FFTs Using Window Functions.

Just keep in mind that if you you increase side-lobe attenuation then you also make your main-lobe wider (worse frequency resolution). Hanning is a good start.

For more info you can refer to literature:

  Harris F.J. - On the Use of Windows for Harmonic Analysis with the
  Discrete Fourier
  Transform
  Heinzel G. - Spectrum and spectral density estimation by the DFT

Also I've noticed in your code that you are normalizing by the number of samples and then taking only first half of your spectrum. That's perfectly ok, but you must also multiply amplitude by 2 in order to compensate for removing of half of your energy in components above $\frac{f_s}{2}$

In case you have more queries you can search for questions and answers regarding windowing and leakage - you will find plenty.


Question based on scaling property of dft Can you please tell how use the scaling property to solve this question?? i am new to dsp subject

Some hints that should help you:

From the scaling property, if we stretch the time axis, we compress the frequency axis.

It is also important to know the periodic property of the DFT: If you extend the frequency axis beyond N samples, the DFT will repeat. Consider a repeating DFT spectrum that extends out 3N, and consider that you now have 3N time samples. You stretched your time waveform from N to 3N, so what occupied 0 to N-1 in the DFT before will now be in the first third of the spectrum (roughly 0 to (N-1)/3). To understand what occurs in the rest you need to be aware of the periodic property of the DFT.

Process some simple DFT's with zero inserts on your own this will also be more obvious to you.


How to model state space for complex valued system correctly in SIMULINK (MATLAB)? When trying to use the default state-space model block, if there is a complex number valued in the matrices, there will be an error

To resolve that, firstly I need to look at pseudo reference model of state space on the internet

Then, create my own state-space block as a new subsystem

u is for input
y is for output
x is the state variable as an output (for sensing use, such as full state feedback)

And here is the inside

Now the problem is that,

The IC (Initial Condition) block cannot simulate what lsim
can done in matlab code, there must be an exponential decay expected
result for the unforced solution, but what I've got is that the
initial condition only happened at the exact t=0 and there is
no exponential decay.  How do I solve it?
The Integrator and IC block cannot operate with complex valued
number, the input and must be separated manually into real and imaginary
part or magnitude and phase part and then joined together again. Is there
any better way to operate them with complex valued number in simulink?
The gain block of A, B, C, and D matrix look ugly, I haven't been
able to change it in the shape of rectangle box as my reference
model above. Is there any way to change it?

Finally, if there is any better solution for my problem, such as pre-builtin complex valued state-space model, or better subsystem, I would be happy to know it

You can split up the real and imaginary part of the state into their own seperate states. Namely by defining $x_r=\mathrm{Re}(x)$, $x_i=\mathrm{Im}(x)$, $A_r=\mathrm{Re}(A)$, $A_i=\mathrm{Im}(A)$, $u_r=\mathrm{Re}(u)$, $u_i=\mathrm{Im}(u)$, $B_r=\mathrm{Re}(B)$ and $B_i=\mathrm{Im}(B)$ then the differential equation can also be written as

$$
\dot{x}_r+i\,\dot{x}_i=(A_r+i\,A_i)(x_r+i\,x_i)+(B_r+i\,B_i)(u_r+i\,u_i)
$$

which when split into their real and imaginary part gives

$$
\begin{bmatrix}
\dot{x}_r \\ \dot{x}_i
\end{bmatrix}=
\begin{bmatrix}
A_r & -A_i \\
A_i & A_r
\end{bmatrix}
\begin{bmatrix}
x_r \\ x_i
\end{bmatrix}+
\begin{bmatrix}
B_r & -B_i \\
B_i & B_r
\end{bmatrix}
\begin{bmatrix}
u_r \\ u_i
\end{bmatrix}.
$$

The output can be expressed using

$$
y=C\begin{bmatrix}
I & i\,I
\end{bmatrix}
\begin{bmatrix}
x_r \\ x_i
\end{bmatrix} + D\begin{bmatrix}
I & i\,I
\end{bmatrix}
\begin{bmatrix}
u_r \\ u_i
\end{bmatrix}.
$$


How do i set the band filter such that it has a pass band of 100Hz to 300Mhz? From this derivation to get the differential of a band pass filter: Id like to find out how we set the upper frequency cut-off and the lower frequency cut-off.

$$i(t)=\frac{v_i(t)-v_o(t)}{R}=C\frac{dv_o(t)}{dt}$$

So we get

$$v_i(t)=v_o(t)+RC\frac{dv_o(t)}{dt}\tag{1}$$

For the high-pass filter we have

$$i(t)=\frac{v_o(t)}{R}=C\frac{d(v_i(t)-v_o(t))}{dt}$$
which gives

$$RC\frac{dv_i(t)}{dt}=v_o(t)+RC\frac{dv_o(t)}{dt}\tag{2}$$

Let's call the time constants of the low-pass and high-pass filters $\tau_L=R_LC_L$ and $\tau_H=R_HC_H$, respectively. For the band-pass filter we need a relation between $v_i(t)$ of the low-pass filter and $v_o(t)$ of the high-pass filter. If we use $v_{oL}(t)$ to denote the output of the low-pass filter, which equals the input of the high-pass filter, we get from (1) by taking the derivative

$$\frac{dv_i(t)}{dt}=\frac{dv_{oL}(t)}{dt}+\tau_L\frac{d^2v_{oL}(t)}{dt^2}\tag{3}$$

From (2) we have

$$\tau_H\frac{dv_{oL}(t)}{dt}=v_o(t)+\tau_H\frac{dv_o(t)}{dt}\tag{4}$$

and (by taking the derivative)

$$\tau_H\frac{d^2v_{oL}(t)}{dt^2}=\frac{dv_o(t)}{dt}+\tau_H\frac{d^2v_o(t)}{dt^2}\tag{5}$$

Plugging (4) and (5) into (3) we finally get for the band-pass filter

$$\tau_H\frac{dv_i(t)}{dt}=v_o(t)+(\tau_L+\tau_H)\frac{dv_o(t)}{dt}+\tau_L\tau_H\frac{d^2v_o(t)}{dt^2}$$

That's much easier in the frequency domain. You simply get
$$V_o(\omega) = V_i(\omega)\cdot \frac{1}{1+j\cdot \omega\cdot R_1 \cdot C_1}\cdot \frac{j\cdot \omega\cdot R_2 \cdot C_2}{1+j\cdot \omega\cdot R_2 \cdot C_2}$$

The first term is the lowpass and the second term is the highpass. To set the frequencies you need
$$R_1 \cdot C_1 = \frac{1}{2\cdot \pi \cdot 300MHz}$$
and
$$R_2 \cdot C_2 = \frac{1}{2\cdot \pi \cdot 100MHz}$$


What does it mean to train and test a data for feature extraction? I know how to do a feature extraction of a data using specific feature like entropy, and energy, using its train and test file to get the accuracy of it. But, i'm still quite confused on the concept of what is the train and test data.

Does train data mean it changes the data after it goes to the feature like entropy?
Then afterwards the accuracy will be checked using the test data?

This is a concept in supervised machine learning.

Train data: Used to train your supervised ML model. This data contains both the input and the desired output, which is compared with the output from your model. Your model is hence trained to reduce the prediction error or cost function.

Test data: The data for which you want to predict the outcome. The desired output is unknown. Based on the training data, your model predicts outcome for the new test data.

You should become more familiar with the concepts of machine learning first, before asking questions here. There is a plethora of resources on machine learning.


Deriving statistics of band limited Random Noise The question is:

Consider a continuous random number with a Gaussian distribution of mean $\mu$ and variance $\sigma$ . The RV is measured from time $t=-\infty$ to $t=\infty$. This time domain signal $x(t)$ is passed through a first order low pass filter with cutoff at $f_o$. What would the statistics of the output be? Would it still be Gaussian? If yes, then what is the new mean and variance?

I am unable to understand where to start with this. Any help is greatly appreciated.

at first read a book on stochastic process or stochastic signal processing.
it's better that we say we have a random process which it's first order distribution (PDF of the process at each instant of time) is independent of time and is a Gaussian with known mean and variance. Now we want to obtain the first order properties of filtered process.

To begin we write the output of filter as a convolution. The output process at each instant of time is a linear combination (weighted average) of input process at different instant of time which all has the same Gaussian distribution, also we know linear combination of variables with Gaussian PDF leads to another variable with Gaussian distribution. So the first order distribution of output must be a Gaussian.

Now taking expectation from both side of convolution, we see expectation of filtered process at each instant of time becomes convolution of filter response and the mean of input process which is a constant over time, so it becomes the multiplication of input mean and the DC response of filter.
To obtain the variance of process around its mean, we have to find the expectation of squared difference of the process and its mean or expectation of squared value of output minus the square of it's mean. To relate the squared value of output to the input we could convolve the filterd signal with its reversed time version at zero lag. At the end we see new variance is the variance of input times square of DC response of filter.


Which Noise Reduction Algorithms Are Used in Commercial RAW Image Processors? I'm trying to guess what noise reduction algorithms are used in commercial processors for raw images from digital cameras. I find this fairly easy to do for the sharpening algorithms (most use unsharp mask as at least one option, for example, and this is clear from the behavior and adjustable parameters), but I know next to nothing about noise reduction. I would like to

identify as well as possible the basic algorithms being used (of course there will be some proprietary optimizations), and
identify the corresponding academic papers that derive them, to learn their inner workings and use as a springboard to the contemporary literature (via tracking citations).

The parameters for the noise reduction algorithms in two common processors, Capture One and Lightroom, are shown below. Here are some further clues.

They're fast. Each works in about a second or less on my modest laptop on 24 MP raw images.
They're old. Lightroom was released in 2007, and I don't think the core options have changed since then.
They operate separately on luminance and color noise.
I don't think they're wavelet-based, since (as far as I know) these applications do not use wavelets anywhere else.

Question 1: Is there a well-known, industry-standard algorithm that meets this description?

Among all raw processors available, the best noise reduction seems to come from DxO's Photolab. (This is of course debatable, but it consistently places at or near the top of comparisons I've seen.) I'd like to try to understand how it operates. A screenshot of its parameters is attached.

DxO describes it as a local algorithm that works pixel-by-pixel:

  “the denoising algorithms (in PRIME) analyze the structure of Raw
  images in depth: more than a thousand neighboring pixels are surveyed
  for each pixel. This extensive exploration identifies similar data
  (for) use (in reconstructing) image information.”

Besides that, we know a few more things:

It's relatively new (released 2014).
It's slow. Users report it takes a minute or more per photo.
It's recommended for very noisy, high ISO photos. DxO and users report no advantage over faster noise reduction techniques for photos with mild noise.

Question 2: What algorithm could DxO plausibly be using?

I realize this is a somewhat unusual question, so thanks in advance for letting me access your expertise.

Common Approaches for Commercial Denoisers
Commercial denoisers are different than what you'd see on most papers. While on papers the results are mostly using objective metrics (PSNR / SSIM) and are evaluated vs. Additive White Gaussian Noise (AWGN) with high level of noise real world images are mostly with moderate level of noise with Mixed Poisson Gaussian Distribution with Spatial Correlation. They are also optimized to look good to the human eyes and not the objective metrics.
Most denoisers out there use 2 main approaches which separates them from what you see on denoising papers:

Multi Scale Approach
The multi scale approach (Commonly done using some Wavelets decomposition but can be done using many non common LPF / HPF Filters or the classic Pyramid Decomposition [Which is a specific case for Wavelets]) is the common way to deal with Spatial Correlation which is common in real world images. The assumption is that for relatively narrow band of the spectrum the noise is white or approximately white which then allows using common denoising methods as seen in papers.

Noise Estimation
Calibration of the algorithm is the key in commercials denoisers. Most of the estimate the noise using single parameter (Noise STD) but the good ones estimate noise curve (Per Scale) which means the noise std vs. luminosity level. It also assists dealing with non standard distributions of the noise. This step is the key for subjectively good results. You can read on it in Darktable Blog Post - Profiling Sensor and Photon Noise. This is usually used with Variance Stabilizing Transform (VST) to deal with non Gaussian Distributions.

Separating Color from Luminosity
The human eyes sensitivity to noise is different in the Colors vs. Luminosity. We're very sensitive to colored nose and less to "Gray" noise. So it is better do the denoising process on Color Channels separately then Luminosity Channel.

Classic denoisier to see which uses the above approaches would be Neat Image.
Common methods used are based on:

Thresholding of a Transformation
Most common transformations would be Wavelets or DCT / DFT based.
Partial Differential Equation Based
Less common and popular but sometimes effective methods. Their evolution started with the Anisotropic Diffusion where later iterations were based on method with Structural / Directional and Steering Vectors.
Non Local Means
Some variant of the Non Local Means denoising method.

Advanced method is based on the merge of the 2 as can be seen in the Block Matching 3D Transform Denoising (BM3D).
DXO Labs Denoisers
DXO Labs historically has used algorithms developed by Jean Michel Morel's Labs and Students.
First iterations were based on the almost vanilla non Local Means algorithms. Later more advanced tooling were added (Multi Scale approach, Noise Estimation, etc..). More advanced iterations are using the concept of Non Local Bayes denoising.
The current DXO Prime algorithm is based on the ideas in the paper - Secrets of Image Denoising Cuisine.
Topaz Labs Denoisers
Topaz Labs with its classic (Non AI) Topaz Denoise has been using non Local Means. This could be concluded from its very specific "Cleaning Artifacts" which are classic to Non Local Means based methods. It had very effective Noise Profiling to performs very well using this method.
Its AI iteration is based on one of the CNN's Denoising methods as the algorithm can not be tuned to a specific dimensions of the image.
Remark
I hope this question will become heavy in details and information from many users of the community.


Camera and laser scanner calibration final approach questions A while ago I posted this question about camera and laser scanner calibration. I've been away from this project for a while and now I need to come back and get a final approach  to calibrate properly this system.

So having Cedron Dawg's as a good answer to get the laser plane and also having the method described in this article, I have the next approach (assuming I got already the camera intrinsic parameters and the distortion is corrected ):

Get the camera extrinsics (R|t) placing a chessboard in front of the
camera on the scanning area.
Get the next equations from this article:

Having the above functions obviously we need a third equation to match the number of equations and variables, so I though to use this answer approach to get the laser plane equation and with this I will be able to solve the three equation system for the world coordinates (X,Y,Z).

So assuming that for every camera frame I've got the laser pixel input image coordinates (x,y) I will be able to transform them to world coordinates (X,Y,Z) with the above equations.

is all of this correct? is there in this approach any mistake?

EDIT

I edit in order to clarify more what I'm trying to do. The next picture illustrates an example about what I'm trying to do:

The object will change in width (B) and height (A) uniformly ( assume laser,camera and target are stationary), so applying the laser I need to measure (for each laser point) height and width changes. So the aim is, for each camera frame, draw a calibrated laser profile of the object.

What would it be the best way to solve this?

Thanks in advance.

The solution I gave in the other answer assumed that the plane in question was fixed relative to the camera location.  Thus a mapping from the 2D pixel locations to a 2D location on the plane could be made.  Using the same calibration points, a reverse mapping can also be made.  These mappings contain the distortion of the lens, the location of the camera relative to the plane, and the perspective effects.  The mappings can be made more accurate by using higher order equations and more calibration points.

Once you have the 2D coordinates on the plane, call them (x,y), you can convert into real world 3D coordinates using a vector parameterization of the plane:

$$ \vec r =  x \vec a + y \vec b + \vec c  $$

Where $\vec a$ is the 3D real world unit vector in the x direction on the plane, $\vec b$ is the 3D real world unit vector in the y direction on the plane, and $\vec c$ is the 3D real world unit vector of the origin on the plane.

If the situation is more complicated than that, e.g. the laser plane moves, or the camera moves, then a different solution is called for.

I have Python code for a generalized mapping solution.  You can contact me at the email address on my profile page.

Hope this helps.

Ced


What does G(1) = 1 say about a system? This is a line from a paper I've been reading:

  The static gain of the closed loop system must be $1$ ($G(1) = 1$) [...]

First of all: I know what gain is, but isn't gain dependent upon frequency? And why does $G(1)$ do the trick of determining the static gain? And what is "static" gain, anyway?

We are talking about a discrete-time closed loop system and $G(z)$ is its transfer function.

Static gain refers to the DC gain. Namely, it would be the ratio of the output and the input under steady state condition.

Due to DC corresponding to $\omega=0$, in the $z$-domain DC would correspond to $z=1$ because $z=re^{j\omega}$, with $r=1$ (i.e. you have to evaluate your transfer function $G(z)$ in the unit circle to get the frequency response).

EDIT:

The frequency response of a system is represented in the Fourier domain, $G(\omega)$. The DTFT of a given discrete sequence can be calculated as:

$$G_F(\omega)=\sum_{n=-\infty}^{\infty}g[n]e^{-j\omega n}$$

On the other hand, the $z$-transform is defined as:

$$G_Z(z)=\sum_{n=-\infty}^{\infty}g[n]z^{-n}$$

It's easy to see that if we evaluate $z=e^{j\omega}$ (the unit circle), then the $z$-transform returns the frequency response.