Is blending of colour lights a linear system? Can we consider recombination process of different light colors a linear system/phenomenon? Like, when you shed two different light sources for example, a blue light and a yellow light on a surface or recombine them through a prism , the resultant light would be green which has a different frequency that the inputs. Producing a new frequency component other than what were in input (Is this assumption correct?) This is contrary to definition of linear systems I assume. So is this process linear, if not, how is it non-linear, I mean the whats the physics?
I've wondered this question from time to time and I know I probably have some wrong assumption some where. I would appreciate if someone clarify things for me or point out my wrong supposition(s).

Your mistake is that you are describing subtractive color mixing, not applicable in this scenario.

If you shed two lights onto a diffuse surface that does not absorb any light, and collect all the light diffusely reflected, the sum will be an additive mixture of the two. Incoherent (non-laser) light can be thought of as noise that has some electromagnetic spectral power distribution. From the signal processing point of view, long-time average power spectra of independent noise signals mix additively when they are summed


When lowpass filtering with a Gaussian kernel, how can $\sigma$ be related to cycles/pixel? I wish to low pass filter an image with a Gaussian kernel. However, I also wish to select a filter that has the property of producing images that have have no greater than N cycles/pixel. Let's say I wish N to be close to, or exactly .025.

By what method is $\sigma$ selected in order to obtain such an image?

The frequency response of a Gaussian is also a Gaussian, see here. The cutoff frequency $\sigma_f$ is defined as the standard deviation of the frequency response. This is given by $$ \sigma_f = \frac{1}{2\pi\sigma} $$ Given that $N$ is the wavelength of the frequency components you wish to remove, you can calculate the corresponding frequency and work out how much you want to reduce it by, work out the corresponding cut-off frequency and filter standard deviation.

A more back of the envelope approach is to say that Gaussians of standard deviation, $\sigma$, remove or strongly reduce components with lateral size $<\sigma$.


scipy.signal.wiener for audio processing Does some have an example of what a Wiener filter (that can be used directly with scipy.signal.wiener) can be useful for, in sound processing (it seems that such adaptive filters can be useful for noise reduction, etc.) ?

I tried with various (noisy + sinusoids) soundfiles (read into an array x) with the command

y = wiener(x)

But it doesn't give such great result... (Normal : I cannot expect something magical with just a simple function that doesn't take any parameter!)

I just wondered what this function could be useful for (sound processing).

This page has some relevant information.

  The Wiener filter is a simple deblurring filter for denoising images. This is not the Wiener filter commonly described in image reconstruction problems but instead it is a simple, local-mean filter. Let $x$ be the input signal, then the output is
  $$
y = \left\{  \begin{array}{cl} \frac{\sigma^2}{\sigma_x^2} m_x + ( 1 - \frac{\sigma^2}{\sigma_x^2}) x & \sigma_x^2 \ge \sigma^2\\
m_x & \sigma_x^2 \lt \sigma^2
\end{array} \right.
$$
  where $m_{x}$ is the local estimate of the mean and $\sigma_{x}^{2}$ is the local estimate of the variance. The window for these estimates is an optional input parameter (default is $3\times3$ ). The parameter $\sigma^{2}$ is a threshold noise parameter. If $\sigma$ is not given then it is estimated as the average of the local variances.


Signal Processing using Fourier Transform How can I derive the fourier transform of

                                   g1(t) = {A cos(2πfc t),−T/2 < t <T/2
                                            0,               elsewhere

and;

                                  g2(t) = cos(t), −∞ < t < ∞

my answer for g1 is;

∫ g(t)exp(-j*2*pi*f*t)dt = ∫ Acos(2*pi*fc*t)exp(-j*2*pi*f*t)dt

   =∫ 1/2A (exp(j*2*pi*fc*t)+exp(-j*2*pi*fc*t))exp(-j*2*pi*f*t)dt

   =∫1/2A(exp(j*2*pi*fc*t)exp(-j*2*pi*fc*t)+∫1/2A(exp        (-j*2*pi*fc*t)exp(-j*2*pi*fc*t))

   =1/2A∫exp(-j*2*pi(f-fc)t) + 1/2A ∫exp(-j*2*pi(f+fc)t)dt

   =1/2A(f-fc)+1/2A(f+fc)

Hey it takes a lot of time in writing equations also I don't know to type in LATEX.so I am sending you image file for second one


DSP using audio IO from PC using C++ This is my first question on dsp. I am looking for something very simple here and I am sure there are libraries that will do pitch recognition and FFT but my requirement is very simple.

I have a PC and I am using the audio input for analog signal feedback and output for controlling the device. I am using windows as the OS.

I want to ensure that I have full control on the audio output and no other program or the OS will stream sounds into the same channel.
I want to output precisely at a known time interval
I want to check if the feedback matches with what is expected of the device.

Any snippet of code or an article to do the same will be useful. I can provide more details if there is something missing in my question.

The best performance you will be able to squeeze out of a PC is with an audio interface that supports ASIO.

To be able to get near 10ms of total delay (input and output) you will also need a very fast machine. Not only in terms of CPU but more importantly in terms of memory performance.

On top of this delay, you are going to have to add the delay in processing your signal and generating the output signal.

1ms (total delay (?) ) on the hardware you are describing is unrealistic.

The way ASIO works will handle issues #1,#2. I am still not really sure what you mean by #3. If you want to use the PC to drive a device with some signal and collect the response of the device to that signal then that is relatively easily achievable with an ASIO enabled audio interface too.

Hope this helps.


Pre-emphasis / De-emphasis filter design I've designed a simple FM communication system and found that high-frequency noise is a real problem.  It appears that this is a well known issue, which has been solved through pre-emphasis of high frequencies at the transmitter and de-emphasis at the receiver.

The standard pre-emphasis filter appears to have the following form:

$$H_{pe}(f) = 1 + j2\pi f\tau_x, \quad \text{where } \tau = \rm75E-6$$

The de-emphasis filter is simply the inverse of that, or

$$H_{de}(f) = \frac{1}{1 + j 2\pi f\tau_x}.$$

I'm afraid I'm not much of an IIR filter designer, but from my research it appears to me that the de-emphasis filter is a relatively simple integrator, and that I can use the Impulse-Invariant method. The pre-emphasis filter, on the other hand, looks to be a bit tougher.  Should I use a bilinear transform approach?  When I tried to use the bilinear command in MATLAB it complained, presumably because there were no poles.

Any tips or hints would be greatly appreciated.

EDIT:
The MATLAB command that I used was:

[z, p, k] = bilinear([j*2*pi*75e-6], [], 1, 20e3);

MATLAB (using version R2011b) returned:

Error using bilinear (line 100)
First two arguments must have the same orientation.

You really can use any filter design scheme that you want. I can't think of any reason why you would need to use an IIR approach, exclusively, so you have a number of options. For FIR filters, you have choices like the frequency sampling and least squares methods (look at fir2 and firls in MATLAB).

For an IIR implementation, you should be able to use the bilinear transform straightforwardly; since the filters you referenced are simple single-pole and single-zero ems, it's probably just as easy to work the math out by hand. One thing to remember is that you'll want to prewarp the frequency axis first to ensure that you get a good approximation across the passband of interest. It might help for you to provide any error that you were getting from MATLAB; the lack of poles shouldn't be problem.

Edit: A few things to address your comments:

I wouldn't say that firls, being a least-squares-based approach, is a step-wise approximation of a smooth curve. Yes, you provide a number of samples of the desired response to the least-squares algorithm, but you should still get a smooth curve out.
I started to work out the math on applying the bilinear transform to the pre-emphasis filter and came across a problem. Its frequenecy response is:
$$
H_{pe}(f) = 1 + j2\pi f \tau_x
$$
Noting that for the frequency response $H_{pe}(f)$, $s = j\omega = j 2 \pi f$:
$$
H_{pe}(s) = 1 + s\tau_x
$$
The bilinear transform maps the $j\omega$ axis in the $s$-plane to the unit circle of the discrete-time approximation's $z$-plane. So, we should be concerned about what happens to the transfer function as $s$ is allowed to vary along the entire imaginary axis. Note that:
$$
\lim_{s \to \infty} H_{pe}(s) = \infty
$$
This means that, for the discrete-time approximation generated by the bilinear transform, as the angular frequency $\omega \to \pm \pi$, the discrete-time filter's magnitude response will blow up to $\infty$, which is probably not what you want.

This problem would present itself for any analog system whose transfer function has a larger-degree numerator than denominator. So it's probably justified why MATLAB issues the following error:

>> [num, den] = bilinear([75e-6 1], 1, 100e3);
Error using bilinear (line 69)
Numerator cannot be higher order than denominator.

(Note that in your example, you passed an empty vector for the denominator; instead, you should pass 1, since the denominator of the transfer function is unity)

So as it turns out, I'm not aware of any analog-mapping techniques that will definitely work for your application. I don't believe the impulse invariant method will work since the pre-emphasis filter isn't bandlimited. You could use one of the FIR design methods, an IIR least-squares technique (e.g. Berchin's frequency-domain least-squares (FDLS) algorithm), or perhaps look into the matched Z-transform method, which I'm not familiar with.


what's the pass band ripple and stop band attenuation of a digital filter? Hi i'm a beginner in signal processing i want to know what'sthe pass band ripple and stop band attenuation of a  digital filter ?
Thanks.

I hope the plot below helps answer your question. Typically I have seen the "passband ripple" and "stopband attenuation" expressed in dB as shown in the picture translating the magnitude of the ripples to dB using $20log_{10}$ as shown.  So the passband ripple is the amount of variation in the amplitude, within the designated passband of the filter, and stop band attenuation is the minimum attenuation level with the designated rejection band of the filter.

Here is a design example showing proper use of the ripple and rejection, along with common techniques used to get a first estimate of the number of taps (in an FIR) that will be needed to achieve the desired specifications.  These estimators have been detailed in other posts under the topic of "How many taps do I need...".


How is downsampling and low-pass filtering done in MS-SSIM? I am just studying the Multiscale Strutural Similarity measure for image quality and going through the original paper [1].

I think I understood the basic ides and "regular" structural similarity, but could someone clarify on how downsampling is done and which low-pass filter is applied when computing the MS-SSIM?

[1] Wang, Zhou, Eero P. Simoncelli, and Alan C. Bovik. "Multiscale structural similarity for image quality assessment." The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003. Vol. 2. Ieee, 2003. Available online.

I assume the question is specifically the decimate by 2 as in this figure:

Decimation is the combination of a low pass filter with a downsampler. The downsampling by 2 designated by the block that has a 2 with a down-arrow is simply done by removing every other sample. This would fold in frequency everything in the upper part of the spectrum so an anti-alias filter is required prior to the down-sampling. This is no different from the requirement for an anti-alias filter prior to an A/D converter. So the ideal filter L would be a low pass filter that would pass everything that would be in the first Nyquist zone defined by the new sampling rate after the downsamping, while completely rejecting all the high frequency components that would otherwise fold in. (Such an ideal filter cannot be realized, but we approach this trading complexity with allowable distortion).

This graphic I have below demonstrates this with a decimate by 4 for a real signal with the first Nyquist zone extending from $0$ to $F_s/2$, where $F_s$ is the sampling rate (and $-F_s/2$ to $+F_s/2$ for complex signals).

With digital sampling (or resampling as in this case) all the spectrum centered on $mF_s$ is aliased to be centered on $0$ in the first Nyquist zone where $m$ is any integer, so here in this example the new sampling rate is $1/4$ of what the original sampling rate, thus we see all the images that are at multiples of this new sampling rate that were part of the original digital spectrum that was at the higher rate.

In this case a multi-band rejection filter is a good solution to minimize resources. Alternate approaches that are also quite efficient are Cascade-Integrator-Comb structures (CIC) and polyphase filters which are detailed in other posts on this site.


How to detect all the connected neighboring pixels which have same color value in opencv I should be able to input a pixel position and get all the same coloured(in my case it should be black) pixels which are connected to it. How to do this in opencv with c++. Simply the output pixels should be connected to each other with color black. FindContours() method does not work as it cannot be feeded with a pixel.

I think you're looking for floodFill. From the documentation:

The functions floodFill fill a connected component starting from the seed point with the specified color. The connectivity is determined by the color/brightness closeness of the neighbor pixels. The pixel at (x,y) is considered to belong to the repainted domain if: [...]
Use these functions to either mark a connected component with the specified color in-place, or build a mask and then extract the contour, or copy the region to another image, and so on.

You're probably interested in this flag parameter:

FLOODFILL_MASK_ONLY If set, the function does not change the image ( newVal is ignored), and only fills the mask with the value specified in bits 8-16 of flags as described above. This option only make sense in function variants that have the mask parameter.


Findin the source direction using cross-correlation function I am new to DSP and I'm working on the following problem.

We have a sound source in the far ﬁeld and a microphone array consisting of two microphones with distance d. We want to determine the source direction by using the cross-correlation of the microphone signals x and y.

I was told that in order to find the angle, we first need to calculte the time delay (in samples) between the microphone signals x and y. I know we can do this by finding the maximum of cross-correlation function R_XY. However, I don't understand why the global maximum of R_XY happenes at this time delay, namely k0.

Let's say we find this time delay. Then how do we calculate the angle? (assuming we know the distance between these two microphones, sound-propagation velocity and the sampling frequency.

With just two microphones you cannot distinguish a unique location on a plane, but you can resolve two possible locations (but your accuracy may be significantly degraded due to reflections/ multipath). This is demonstrated in the graphic below:

Given one microphone and one resolved range (from time delay and sound propagation velocity), the possible location of the source is anywhere on a circle (on a 2D plane, or anywhere on a sphere in 3D space).  With two microphones, this can be reduced to two points (on a 2D plane, or a circle in 3D space where the two spheres overlap).

To uniquely determine the location of the source on a 2D plane, 3 microphones would be required. But if you have further information, such as you know the source is in front of you such that you only need to resolve azimuth angle from -90° to +90°, then one of the two points that are resolved can be neglected.

Therefore the solution is solving the for the x,y location given the simultaneous equations for two circles on a 2D plane:

$$(x-x_1)^2 + (y-y_1)^2 = r_1^2$$
$$(x-x_2)^2 + (y-y_2)^2 = r_2^2$$

As a reminder if needed such a solution is worked out here: http://www.analyzemath.com/CircleEq/circle_intersection.html

By the way, this is exactly how location is determined with GPS (triangulation) as depicted in the graphic below showing how a unique location on a plane can be determined with known range to three satellites:

The graphic says "Simplified" as it is location on a 2D plane, but the explanation above using spheres instead of circles applies to determine location in a 3D space. Give the knowledge that the satellites are above, only 3 satellites would actually be required if we had accurate time as the other location in 3D space is in outer space so can be neglected. However, since we do not have sufficient time accuracy in the receiver, a 4th satellite is needed to determine that final unknown (thus with 4 satellites and 4 ranges to each we end up with 4 equations with 4 unknowns which can all be solved: x, y, z location in 3D space and time t).


VHDL fixed point implementation of IIR filter I have a filter described by the following transfer function H(z):

$$H(z) = \frac{156 - 156z^{-2}}{16384 - 32443z^{-1}+16073z^{-2}}$$

How will this be implemented with the first term in the polynomium being different from 1?

Normally i would implement the filter with something like:
$$y[0]=156\cdot x[0]-156\cdot x[2] +32443\cdot y[1]-16073\cdot y[2]$$

but in this case, it is missing the term:
$$-16384\cdot y[0]$$

How may this term be inserted? I have tried just inserting it, inserting and shift the y's by 1. I just get noise doing this (looks like overflow).

It is to be implemented in VHDL on a FPGA.

Filter output has to be divided by the first term of denominator.
In this case, division by 16384 can be effectively done by logical shift 14 bits to the right.


Cramer-Rao Lower Bound In estimation problems, we may use Cramer-Rao Lower Bound (CRLB) to evaluate the best performance. But if there is no unbiased estimator can attain CRLB, what is the meaning of CRLB?

To clarify the question, I present a simple question here:

$\alpha$ is selected as the evaluating indicator of system $A$. We can obtain the CRLB of $\alpha$, but there is no estimation $\hat \alpha$ which can reach CRLB. Is it reasonable to use the CRLB of $\alpha$ to evaluate the performance of system $A$.

CRLB is usually used to assess the performance for an unbiased estimator because it is seen as the best variance you can achieve. If I understand your question correctly, you are asking if it is fair to use CRLB to assess the performance for a biased estimator? Well, this may not make much sense because you may come across situations where the biased estimator has better variance than the CRLB!

For example, Suppose the true value = 0, you can see one estimate is unbiased with large variance but the other has a bias with smaller variance.


Question about correlation of two signals I'm working on some practice examples and in one of them I'm asked to find the convolution and the correlation of some signal pairs given. I know that correlation , given two signals x(n),g(n), is given by the following :
$$Φ_{xg}=x(n)\ast g(-n)$$ where asterisk is convolution. Now is this the same as $$g(n)\ast x(-n)$$

I think I read this somewhere $$Φ_{xy}(n)=Φ_{yx}(-n)$$ So the 2 convolutions give different results. Shouldn't the question specify which correlation I am to calculate? $$Φ_{xy}\\Φ_{yx}$$

You are correct in noting that cross-correlation is not commutative. If the question says "correlation between two signals $x$ and $g$" then it is usually assumed that you go in that sequence and compute $\Phi_{xg}$.

Or you can choose to say that the question and ambiguous compute both $\Phi_{xg}$ and $\Phi_{gx}.$


Specs for creating sample set? I would like to create a set of audio samples for playback as a midi-controlled instrument. Thus: a sample for every key of the piano range, with multiple volume timbres for each note. I am able to create recordings at various specifications. My problem is finding a specification as to how a sample set is packaged so that it can be used by a standard DAW. Searching so far has not turned up any leads.

If there is a appropriate in another forum to try, I'd be happy to learn of it.

SoundFont 2.0 (SF2) sampled sound banks can be used in any digital audio workstation (DAW) software either natively or via a plug-in virtual instrument. It is an old format from the 16-bit Sound Blaster days but still widely supported, and satisfies your requirements of having key and velocity ranges. More recently there is the SFZ file format that is manually editable as a text file. You tell in the SFZ file which WAV files are to be used and how. There are also SF2 and SFZ editors and SF2→SFZ converters available.


Physical meaning of rippled spectrum What do ripples in the spectrum of a bandpass filter physically mean? What do they tell us about the filter? Please check the picture:

The filter should look smooth, like a resonator, such as

Ripples in a spectrum are an indication of echoes in time.  Specifically consider how the Fourier Transform of a sinewave is two impulses in frequency. Similarly, given the reciprocity of the Fourier Transform, a sinewave in frequency (in other words, ripples in the spectrum) is two impulses in the time domain. The Fourier Transform of the frequency response is the impulse response. A system that has an impulse response consisting of two impulses at different locations in time will output the input combined with a delayed copy of the input (echo).
In practice this can be caused by multipath delays for a signal propagating between a transmitter and receiver, reflections from poor matching on a high frequency (microwave) circuit, or other equivalent effects that result in multiple copies of the signal.


First derivative analog filter I'm reading about fault detection via signal processing in time domain. One possibility is to check that first derivative of the signal is in some predefined bounds. The text says that to obtain the first derivative of the output signal $y(t)$, I can use this analog filter:
$$F(p)=\frac{p}{\frac{1}{\tau}p+1}$$
for some sufficiently large $\tau$. Why is that? I don't  know much about analog filters. I only know that Laplace transform of the first derivative of a function $f(t)$ with zero initial conditions is $pF(p)$. I also understand that $X(p)=\frac{1}{\frac{1}{\tau}p+1}$ is Laplace transform of exponential decay function $e^{-\tau t}$. Could you please explain to me why the filter needs to be designed like that and how does it work physically?

Thank you.

My wrong reasoning: If $u(t)$ is the input signal and $y(t)$ is the rate of that signal, then
$$
y(t) = u(t)'. \\
$$
The transfer function would then be
$$
\begin{align}
Y &= pU \\
F &= \frac{Y}{U} = p.
\end{align}
$$

I will answer this with the understanding that $p$ is the Laplace variable (which I will call $s$).

As you stated, differentiation in the time domain corresponds to multiplication by $s$ in the Laplace domain. So it seems that we could get what we want by just choosing $F(s)=s$; then the output is the derivative of the input, i.e.,

$$y(t) = \frac{du}{dt}$$.

So why the $\tau s+1$ in the denominator?

This corresponds to a lowpass filter with bandwidth $1/\tau$. The output of $F(s) = s/(\tau s+1)$ is the derivative of the input signal after it has been smoothed by a lowpass filter with the desired bandwidth (or alternatively, this is the same as taking the derivative of the input signal and then smoothing it with the lowpass filter). This prevents the derivative operator from picking up high-frequency components of the input signal and amplifying noise.

I don't know much about fault detection, but I can give a simple example. Say you're trying to detect when $u'(t)$ exceeds a fixed threshold. If you had no lowpass filter, small-amplitude (but high-frequency) noise in your measurement $u$ could trigger your detector even when it really shouldn't. But by setting $\tau$ you can filter out components of $u$ above a certain frequency ($1/\tau$) so they don't affect your derivative estimate.


Frequent update of the Kalman Filter covariance matrices I have inherited some code that does vehicle localisation based on two main inputs: sensor input from the vehicle odometry and map localisation based on matching camera images to road features. Both tracking routines works reasonably well, and the outputs from both are fed into a Kalman Filter (the OpenCV one is used) with a reasonably good output under normal circumstances. However, as we round corners, we get a lot of lateral drift and lose visual tracking, so to compensate (?), we double or so the covariance for the X, Y and Yaw values for the measurementNoiseCov, as OpenCV calls it, resetting them after we get back on track; depending exactly on how we fudge the values, the match to Ground Truth improves to some degree, but I think we're just stabbing in the dark and getting lucky occasionally.

Does this make sense?

We also further apply a Jacobian matrix to the X and Y covariance values (sorry, not familiar with the markup:

Jac = [cos(yaw) -sin(yaw) sin(yaw) cos(yaw)]
measurementNoiseCov = Jac * measurementNoiseCov * Jac.t()

Does this also make sense? I've tried reading up various papers, but I haven't really been able to make much sense of them.

The term "Kalman Filter"  tends to be used in genetically and there are a number of variations suitable for different applications.

The documentation in OpenCV for the Kalman Filter describes a strictly Linear type which they say can be modified to be an Extended Kalman Filter (EKF). Your Jacobian implies that your problem is nonlinear, and the OpenCV version is exclusively a linearization.  Many problems are solvable using this approach, but it works better, the smaller your time step (in other terms updated more often). You can also inflate your state covariance.

Another approach, where the state transition $\bf{x_k} \; \Rightarrow \;\bf{x_{k+1}}$ is nonlinear, and using an integration technique, typically Euler Integration is used on the State Evolution equation. The nonlinear state covariance update uses the Jacobian linearization. It doesn't look like you can use the OpenCV KF, for this kind of EKF.

One issue with the KF in general is choosing the number of states you use in your filter. A very common problem in applications is when there is a large change in state, like turns, the Filter exhibits a lot of lag.  If an additional state was present to account for the change (acceleration in Kinematic models) the lag would be much lower, but the the higher dimension can make the Filter more sensitive to small disturbances.  In some cases, this isn't a problem and in others there is the Interacting Multiple Model (IMM) or as a friend of mine called it "Dueling Kalman Filters". The OpenCV doesn't look like it provides the innovations, which probably means that IMM is not a candidate with an all linearized EKF.

The essential problem when you have a jump and don't account for the possibility of a jump, is that the Filter places more trust on the predicted state and the change looks like noise.  It wants (as an anthropomorphic analogy) to keep the state level and steady.  In many cases, the error accumulates slow enough for the Filter to adjust to the jump.  In other cases, track is lost. State Covariance inflation is one way to have the Filter weigh the new measurements over the predictions based on the old state.

One way to avoid linearization problems and approximate integration is to used the Unscented KF, which is often a good solution. The state dimension problem can creep in with this filter a well.

The sensitivity to jump lag manifests itself differently to different forms of the KF and update rate dependency is present in all forms.  Kalman Filters need a lot of tuning and testing in any case.

Step size or update rate implicitly effects time integration error. Small steps, less error. Nonlinear integration, less error.


Extracting a Varying DC from a Signal I have a signal that has both AC and DC components and both a re varying with time. How can I extract the varying DC wave from the signal?

eg: A PPG (Photopletysmography) signal extracted from a healthy patient.

DC, by definition, is the non-varying part of the signal. To find the DC component you simply average the entire signal.

  How can I extract the varying DC wave from the signal?

This sounds more like you want to extract the low frequency part of the signal. The most basic filter for this is a moving average filter. However, I would recommend a smoother filter, such as a Gaussian filter. Start with $\sigma$ about 1/6 the period of the low frequency oscillations you are interested in and modify to suit.


Orthogonality of filter impulse response to its even shift I meet this problem but still do not know how to solve it.
Could you guy give me some guides?

Upsampling by 2 ($U_2$) followed by filtering by $g$, with operator $G$

And given: $<g_n,g_{n-2k}>_n = \delta_k$
(Filters with impulse responses orthogonal to their even shifts)

Prove that:
$$
I = U_2^*G^*GU_2
$$

Thank in advance.

Since this is clearly homework, let's help you with hints:

The set of  tricks for such problems is typically pretty limited; try (and that's really what solves the issue)

Multiplying the equation with elements from your right-hand term from left and/or right, so that identity matrices form.
Writing your matrix/operator property down as structure of your matrix
applying transposes, and hermitians until things happen.


Does chirp have constant magnitude frequency response? Pg. 223 claims so, yet my results via DFT differ:

Is the textbook wrong?

My attempted explanations: (code)

DFT vs DTFT: "frequency response" is computed via latter. Still, DFT should resemble a sampled DTFT.

DFT time-domain periodicity, whereas DTFT assumes aperiodic, or "repeats at infinity" with infinite zero-padding.

To address each, I try greater N, and zero-padding - below. Zero-padding appears to correct phase (quadratic if unrolled), and more samples tend to flatten the magnitude for an ever-growing portion of frequencies to the right.
I figure, in the limit N -> inf, the amplitude spike has zero width (like in Gibbs) - but this appears contradicted in the "large N long padding" case, where a nontrivial portion of the amplitude decays with oscillations. Further, the left peak appears to scale with N, behaving more like an impulse in the limit, which won't yield zero energy as in Gibbs phenomenon.

Update: turns out magnitude doesn't spike, but rather decays exponentially, which is far from the expected horizontal line - and, the phase is linear:

Note: see comments below accepted answer for further info.

The book is not wrong, but it does present the concepts on LFM in a clunky manner and can be misleading. The book presents the analytical expression for the LFM spectrum, which is an approximation. It also plays with the plot views and most likely unwraps the phase angles, which is usually required to see the phases you expect.
Usually when you're introducing LFM, you'll show the modulated pulse itself as well as the phase progression in the time-domain. The analytical expressions in the time domain is all you need to observe the linear frequency and therefore quadratic phase nature of LFM. Doing this in the frequency-domain just tends to introduce more confusion. An example of the time-domain LFM pulse and its phase is shown below.

When deriving the expression for the Fourier transform of an LFM pulse, you do indeed yield a magnitude of 1 over the bandwidth of the pulse. This is intuitively satisfying because you have the same contribution from each frequency over the bandwidth.
Confusion does occur however when one goes to plot this if they expect a constant frequency response. With any practical LFM spectrum plot, even with very long pulse widths, one should expect a ripple effect which you have already identified. The quadratic nature of the phase is still captured in the DFT. The spectrum of the LFM pulse from above is shown below.

I haven't tried your code, but it might be that you just need to zoom in the proper areas and unwrap the phases to see what you want. Provided is the MATLAB code to produce the plots to help you in converting it to python.
%% LFM - Time and Frequency Domain

% Sampling
Fs = 50e6;

% Pulse parameters
tau = 50e-6;
bandWidth = 10e6;
alpha = bandWidth/tau;

% Define waveform
t = 0:1/Fs:tau - 1/Fs;
fmcwPulse = exp(1i*pi*alpha.*t.^2); % Complex transmitted LFM waveform

% Plot
figure;
subplot(2, 1, 1);
plot(real(fmcwPulse));
xlabel("Samples");
title("LFM Pulse - Real Part");

subplot(2, 1, 2);
plot(unwrap(angle(fmcwPulse)));
xlabel("Samples");
title("LFM Phase");

figure;
subplot(2, 1, 1);
plot(abs(fftshift(fft(fmcwPulse))));
title("LFM Spectrum Magnitude");

subplot(2, 1, 2);
plot(unwrap(angle(fftshift(fft(fmcwPulse)))));
title("LFM Spectrum Phase");

Update
Modifying the code above so that $\tau = 1 \space s$, which is relatively long, yields a spectrum closer to the ideal flat spectrum that one would expect analytically. The spectrum is shown below.


Measuring the contrast of an image For a project, I want to measure the contrast of a global image.

I want to compare different images to get the image with the best contrast.
I know about the Weber contrast, but I need a normalized contrast value (in the range between 0-1 or other).

Is there a simple method to do this in openCV?

One simple way for quantification of contract that I can think of is through use of image histogram. Following is my suggestion

Compute Histogram of the Image
From the counts compute entropy

If you just want to try it out you can use the matlab inbuilt function http://www.mathworks.ch/ch/help/images/ref/entropy.html

You can use the entropy value of the histogram as a measure of contract.


Simulate BER curves for OFDM with convolutional coding I am tasked with generating performance curves for an OFDM transmitter/receiver, and I am unsure how to calculate the appropriate AWGN power, and where exactly the apply it.

The transmitted packets contains typical non-payload portions, such as a synchronization marker, short training sequence, long training sequence, cyclic prefix, pilot subcarriers, and guard band subcarriers. The models have quite a bit happening in terms of real-world effects: Carrier Frequency Offset simulation, multipath simulation, channel-specific effects simulation.  These are used to exercise the receivers ability to handle these real world effects, using CFO correction as well as Equalization.

Furthermore, the payload bits are applied to a 1/2 convolutional encoder/decoder, so the IQ symbols look like QAM-4 (2-bits per subcarrier), but after the decoding happens, it's an effective 1-bit per subcarrier.

How can I apply noise the transmitted waveform to correctly generate a BER vs EbNo plot?

Look at one OFDM symbol at the point where it is transmitted, i.e. right before it enters the channel. Figure out its average power, $P$, or its average energy, $E = PT$, where $T$ is the symbol's duration. Also, figure out how many information bits, $I$, are carried by the symbol.

In your simulation, $E_b = E/I$.

To get a desired $E_b/N_0$, calculate $\sigma_n^2 = N_0 / 2$. Then, add AWGN noise to the symbol right before it enters the receiver; i.e. after all other channel effects such as multipath. The noise should have variance $\sigma_n^2$ in each of the I and Q branches.


Linear model of the modified error function We can write the error function ($E(w)=1/2\sum_{n=1}^{N}\{y(x_n,w)-t_n\}^2$) as a linear model using its partial derivatives. Is it possible to do the same thing about the modified error function?

The modified error function is (page 10- pattern recognition by Bishop):

$$E(w)=1/2\sum_{n=1}^{N}\{y(x_n,w)-t_n\}^2+\dfrac \lambda2||w||^2$$

Yes, you can do that.

$$E(w)=1/2\sum_{n=1}^{N}\{y(x_n,w)-t_n\}^2+\dfrac \lambda2||w||^2$$

$$\begin{align} \frac{\partial E(w)}{\partial w} &= \frac{\partial \left( 1/2\sum_{n=1}^{N}\{y(x_n,w)-t_n\}^2+\dfrac \lambda2||w||^2 \right)}{\partial w} \\ &= \frac{\partial \left( 1/2(y({\bf x},{\bf w} )-{\bf t})^T(y({\bf x},{\bf w} )-{\bf t})+\dfrac \lambda2{\bf w}^T{\bf w} \right)}{\partial {\bf w}} \\ &= \frac{\partial \left( 1/2(y({\bf x},{\bf w} )-{\bf t})^T(y({\bf x},{\bf w} )-{\bf t})\right)}{\partial {\bf w}}+\frac{\partial \left( \dfrac \lambda2{\bf w}^T{\bf w} \right)}{\partial {\bf w}} \\ &= \left( \frac{\partial y({\bf x},{\bf w} )}{\partial {\bf w} } \right)^T \left( y({\bf x},{\bf w} )-{\bf t} \right) + \lambda{ \bf w }\end{align}$$


How to evaluate a close-loop control during transient? I have implement a close-loop control in order to perform the tracking of the desired signal R.
The basic scenario shall be summarized as follow:

I have to evaluate the performance of my control G with a moving average value of the Error E.
When there is no transient the error E shall be less then a certain threshold.

But, how to relax the threshold during transient ?
Addition info:

Frequency of R is limit to 2 HZ
Rising time is 200ms
Control G has a 1 Hz band (but this is not a problem)
I can acquire R E U and Y
Real-time FFT of R is not possible during execution.

My idea is to find a threshold as function of r signal which is small when reference is stationary (constants from at least 200 ms) and high when the reference is fast changing.

Thank you all
Regards

In order to evaluate the performance of a controller during transients, you can use different metrics. Two metrics you could use are

ITAE= $$\int_0^{\infty} t*|e| \,dt$$

ITSE = $$\int_0^{\infty} t*e^2 \,dt$$

where e is the difference between your reference and the measurement and t is the time since the reference change.

Basically, you want to minimize these metrics. The basic idea is that an error in the beginning won't impact the metric much while an error 1 second after the reference change will impact your metric much more.

There are probably other metrics you could use. You could also use stabilization time at 1% or 5%.

Edit : If your signal is cyclical you could limit the integral upper bound to a limited number of 2-Hz cycles.

ITAE= $$\int_0^{NT} t*|e| \,dt$$

ITSE = $$\int_0^{NT} t*e^2 \,dt$$

where N is a whole number and T is the period of your reference signal.


Per Window Normalization I have PPG signal that should be normalized, with offline processing I'd just do $z$-normalization,
$$
\displaystyle \frac{x- \mu_x}{\sigma_x} \quad \text{where } x \text{ is the signal.}
$$
The problem is that I get the data in chunks and I'm not sure what would be the effect on the signal when normalizing only chunks of the data at any given point. Any ideas or suggestions on how to approach this would be highly appreciated.

Thanks,
Jack

In several application the mean is removed by using a high pass filter. A first order high pass filter has the following structure:
$$y[k] = x[k] - x[k-1] + \beta y[k-1]$$
or
$$y[k] = \beta x[k] - \beta x[k-1] + \beta y[k-1]$$
where $0<\beta<1$ determines the steepness of the high pass filter, larger $\beta$ leads to a steeper/sharper filter.

NOTE: Multiplying $x[k]$ and $x[k-1]$ by $\beta$ (as in the second equation) ensures that the gain per frequency is at most one. Due to the next normalization step, using this $\beta$ for $x$ is redundant.

Normalization of the signals energy can be performed with a recursive estimation of the signal energy:
$$ P_x[k] = \alpha y[k]y[k] + (1-\alpha) P_x[k-1] $$
where $y[n]$ is the zero mean signal and $0<\alpha<1$ is a forgetting factor. Small values for $\alpha$ lead to large time-constants of the recursive filter.

The normalized output output $z[k]$ is calculated as follows:
$$z[k] = \frac{y[k]}{\sqrt{P_x[k]}}$$

The signal $z[k]$ is a zero-mean signal with unit variance. You can scale your variance to a desired level $\sigma^2$ by applying the gain $\sqrt{\sigma^2}$ to signal $z[k]$.

NOTE1: A signal in the rang [0,1] always has a non-zero mean if the signal is not zero all the time.

NOTE2: In order to guarantee the signal levels to be in a specified range such as [-1,1], either a compression technique has to be applied or you have to know the original abolute maxima and minima such that you can scale to the desired range.

NOTE3: The variance in combination with the distribution only tells something about the majority of the samples.


Adding channel effects to a signal I have a a channel matrix "H" that is circulant.
I have data blocks.
I want to add the channel effects to the signal.

When H was only a vector of channel coefficients I would say:

%Going Through The Channel
After_channel= filter(H,1,Data);

but now that H is a matrix the line above wouldn't work.

I'm not sure what to do

When you use the channel matrix (it is probably a convolution matrix), you just need to multiply it by the signal. Note that the output from this operation will not be equal the output from the filter function. The filter function generates an output from the same size of the input signal, so if your signal has a length of N, just the first N samples obtained from the multiplication will be the same. See the code posted by @Harris.


Other time-frequency-plane tiling than STFT, DWT, ConstantQ-Transform: multiresolution STFT? It is known that

a) the STFT gives a rectangular tiling of the time-frequency plane

b) the Wavelet transform gives a non-linear tiling (better frequency resolution for low-frequencies, and better time-domain resolution for higher-frequencies)

c) Constant-Q transform (such as NonStationaryGaborTransform) have a logarithmic scale for frequency bins (instead of linear with STFT) and have a time-frequency tiling like this (y-axis logarithmic):

Question:
Is there an transform like this:

i.e. like a normal STFT, but for which the FFT-size would change for successive time frames.
Example: if there is a transient, the FFT-size is small (512) to keep good time-domain resolution, then a few time-frames later the signal is rather stationary, so the FFT-size is higher (8192) to have good frequency resolution.
It's like an adaptative / multiresolution STFT.
Last but not least: of course it's always possible to perform FFT with different window-length on successive time-frames, so the forward transform I'm describing is probably easy to do.
What I'm asking here is a transform which is invertible, i.e. in which we can compute an inverse transform to recover the initial signal.
Does this exist in Matlab or Python?

Your time-frequency grid is Mondrian-shaped: essentially, rectangles supported on dyadic splits of the time or the frequency axes. Hence, you can easily start for any reversible time-frequency tiling, and further slit any of its rectangle onto another invertible tiling. Some call that hierarchical, or nested time-scale/time-frequency decompositions. As long as each is invertible, the combination remains revertible.

Furthermore, dyadic and power-of-two uniform decompositions can be turned into the other, only by grouping coefficients. This was used for instance in A DCT-based embedded image coder, 1996: an eight (or $2^3$) multiband filter-bank DCT was turned into a 3-level dyadic structure.

So, combining dyadic wavelets and $2^K$ windowed filter-bank, you can get invertible schemes as wished.

One issue though: how does one choose how to nest those different embedded decompositions? Seen as wavelet or local cosine, lapped-transform packets, entropy concepts could be used, but the diversity can be huge.


Creating odd harmonics without the fundamental I'm interested in creating a series of odd only harmonics, without the fundamental, for audio use.

As an example the input could be a sine wave at specific frequency $f$ and output is $3f$, $5f$, $7f$, etc. - but without the original $1f$.

Since the $1f$ (in a real audio system) is unknown and the input is multi-tone, just filtering the $1f$ out is not an option.

Any idea for replacing an input of $1f$ with an output of $3f$, $5f$, $7f$, etc. ?

Edit: Additional information from the deleted post

Example: Let's assume the fundamental to be anything possible between 100 and 4000 Hz. Whats more, what if the input is multitone (as most music is...). 50 and 425 Hz input simulatenously ??

Im looking for a mathematical function which can recreate any fundamentals odd harmonics, in a multitone input...

Im looking for a mathematical function which can recreate any fundamentals odd harmonics, in a multitone input...

(emphasis by me)

What you ask for is not possible. If you demand multiple tones and complex tones as they appear in music, then your problem leads to a conflict. A complex tone comes with overtones. Even if you have a process that adds odd harmonics to each sinusoidal component, then the even overtones of the fundamental frequency will contribute even harmonics. That is because an odd number multiplied with an even number is again even.

So circumvent that, you would have to identify the fundamental frequency of complex tones first, leading you directly to the problem of polyphonic pitch detection which remains mostly unsolved.

Maybe you should explain what you're really trying to accomplish, as in what the application is that you have in mind. With that information it's quite likely that we can come up with a relaxed set of requirements that can actually be implemented.


Confused on the difference between the frequency spectrum of an entire song, and the frequency spectrum of a point in time I am very much a beginner in this field - but find it really interesting. However I am a little confused on a certain area of knowledge.

If I have understood correctly:
At any point in an audiosignal, for instance when someone is strumming a D-Major on a guitar, i can reproduce that exact sound in that moment by looking at the frequency spectrum for that time, do a Fourier Transform, and recieve a periodic signal in the time domain presented as a fourier series. If I play that signal, I will hear the sound that was in that exact moment of a song, just as a never ending note. In this scenario the frequency spectrum tells me what pitches were played- and the frequency domain will change for point in time.

However, If a look at a three-minute-song in the time domain, I dont quite understand what the fourier-transform will tell me. If I have this audiofile, and I view it in the timedomain, what will it tell me? (view picture)
And if I transform it to the freq. domain, what will it tell me?

As far as my logic goes, I wont ge able to tell any of musical features just of an amplitude value in the time domain.

(And the confusion continues when we sample an audiosignal, and all we get is an amplitude value for a given time. How can any musical information be there?)

Thanks for answers.

As far as my logic goes, I wont ge able to tell any of musical features just of an amplitude value in the time domain.

That's wrong. The amplitude-over-time actually is the way the song is played by your sound card (we call that PCM). So, this is your song.

What you really need to understand that your analog audio (i.e. the variations in air pressure over time) is being sampled and stored digitally, and reconstruction happens based on these samples.

Aside from that, you're absolutely right, the Fourier Transform and the time domain signal are equivalent representations of the same signal.

Regarding the question what a graphical representation will tell you: That's quite impossible to answer, as it depends on you ;) Now, there's a few basic things that one can note: For example, if you look at a plot of the time domain signal, and you see high amplitudes in a certain range of time, you can assume that it's loud during that time. Same goes for the absolute of the Fourier Transform: If the amplitude spectrum is large for a range of frequencies, then you can assume that a loud tone of that frequency range is audible.


Can all systems be described by differential/difference equations? All LCCODEs (linear, constant-coefficient ordinary differential equation) with initial rest conditions describe causal LTI systems (is this a correct statement?)

Is the converse true, i.e. can all causal LTI systems be described by LCCODEs with initial rest? If not, what are the necessary/sufficient conditions?

It seems trivial that arbitrary differential/difference equations describe some sort of system (is that also true?)

Again, is the converse true, i.e. can we describe any arbitrary systems with (possibly
nonlinear/variable-coefficient/partial, without restriction on initial
conditions) differential equations? If not, what are the necessary/sufficient conditions?

No the converse is not true in general.

Take for example the discrete-time ideal lowpass filter with impulse response

$$ h[n] = \frac{ \sin( \omega_c n) }{ \pi n } ~~~,~~~-\infty < n <\infty$$

which describes an LTI system but it does not correspond to a difference equation of any kind. Indeed $h[n]$ is derived based on the inverse discrete-time Fourier transform of the filter frequency response;

$$ H(\omega) = \begin{cases} {1 ~~~,~~~ |\omega| < \omega_c \\ 0 ~~~,~~~ \text{ otherwise } }\end{cases} $$

$$h[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} H(\omega) e^{j\omega n} d\omega $$


Conceptual question on FFT/IFFT (IFFT existence) I was reading  "Discrete and continuous Fourier transforms: analysis, applications and fast algorithms" written by E.Chu and, at some point, I found something that I could not completly understand.

  If $x(t)$ is not continuous, then it is possible that 2 different functions share the same transform $X(f)$. Therefore, the function $x(t)$ cannot be uniquely determined by inverse trasforming $X(f)$ unless it is continuous or it satisfies 2 more conditions:

  1- $x(t)$ has only a finite number of maxima and minima on any finite interval

  2- $x(t)$ has on any finite interval at most a finite number of discontinuities, each of which is a jump discontinuity.

  In the latter case, the inverse Fourier transform of $X(f)$ produces $\hat{x}(t)$, which agrees with $x(t)$ at every $t$ at which $x(t)$ is continuous, and $\hat{x}(t)$ equals the average of the left-hand limit $x(t^-_{\alpha})$ and the right-hand limit $x(t^+_{\alpha})$ at every $t_{\alpha}$ at which a jump discontinuity occurs.

So, here my doubts.

First of all, I am not sure I got the meaning of the statements related to the 2 conditions listed. This is probably due to the fact that I am familiar only with signals coming from structural tests (mainly accelerations), therefore I think it could be really helpful if someone of you can provide me some example I can focus on to understand the concept (for example, how can be possible to have infinite maxima and minima in a signal? And what are the discontinuities it is talking about?).

Second, measured signals (for example, the aformentioned acceleration measurements coming from a test) are not continuous (since they are sampled at a given sampling frequency), so if I want to apply the IFFT algorithm to retrieve the time signal does it mean I have always to check if the 2 other conditions are satisfied? And how can I check them? Visual inspection? Or do I have to look for some specific features the signal must have?

Are there classes of signals which always satisfy the 2 conditions?

The conditions you mention only apply to the continuous-time Fourier transform, and they only have theoretical relevance. All signals occurring in engineering applications will satisfy these conditions. But if you want an example of a function with an infinite number of maxima and minima on a finite interval, here it is:

$$x(t)=\sin(1/t)$$

This function has an infinite number of minima and maxima on the interval $(0,a]$ for any real-valued $a>0$.

And as for your measured signals, they are all discrete-time signals of finite length, and the above conditions do not apply. So you don't need to worry about them when using the DFT/FFT, because the DFT is defined by a finite sum, unlike the continuous-time Fourier transform, which is defined by an improper integral.


filter design steps and process I started to learn MATLAB and reading paper.i have faced difficulty where to start, how to start in the field of biomedical signal filtering (fNIR)?

Your question is overly broad and fundamentally impossible to answer because one needs preparation and no one knows how prepared you are.

Matlab and signal processing are 2 different things. It is used a lot by people who do signal processing but they are not the same.  Learning Matlab if you have some programming background and some linear algebra is not hard.  I’ve seen 60 year old die hard FORTRAN programming engineers pick it up in 2 weeks.  They understood the numerics and it was just a matter of reexpressing their knowledge.

Using a tool you don’t know how to use on a problem you don’t understand is very difficult. Did I do it right? Is it the code? Is it the algorithm?  Often it will be both.

FDTOOL is considered easy to use.

The Understanding Digital Signal Processing by Lyons is recommended.

Look through the old posts on this sight, particularly those that are “frequent”.

Some people learn in a line. Other people learn in a circle. They revisit concepts after learning new ones. Circles have no start and end. I think everyone is a bit of both.


non-consistent normalization problem in fft I have a csv file of data sampled with Ts=1ns which looks like this:

This signal is a step response of some system which responds to a step of value 1.
I'm trying to get the impedance profile of the signal by using
$$
Z\left( f \right)=\frac{V\left( f \right)}{I\left( f \right)}=\mathcal{F}\left\{ \frac{1}{A}\frac{\partial }{\partial t}v\left( t \right) \right\}
$$
The signal is noisy and needs filtering, so i'm separating the droop part from the non-droop part and filtering them independently using savgol filter.

Than, i'm applying FFT on each of the signal parts above using the following code:
n = 1024
from scipy.fftpack import diff
frqs = np.fft.fftfreq(n, d=1.0 * 1.e-9)[range(n // 2)]
Y = np.fft.fft(diff(data_y, order=1), n=n, norm='ortho')  # fft computing and normalization
Y = abs(Y[range(len(frqs) // 2)])
return frqs, Y

Next, I take the high frequencies from the FFT applied on the droop part and the rest from the non-droop part and I stitch the results and plot them in log scale.
We get a result with expected shape as below:

Now, the problem arise when we change the droop part to be longer (and thus, the non-droop part to be shorter) as seen here:

The output of this looks like this:

The noisy shape is expected, but the amplitude is changing drastically and I can't figure out why.
can someone point what the problems are are help with a solution?
Thanks!

I believe you need to scale by the sample period $\Delta T$ properly in your differentiator (divide by $\Delta T$) if your time units are in seconds not samples.

Your differentiator can be arrived at by using the backward Euler approximation of a continuous time differentiation.

The mapping from s to z using the backward Euler

$$s \leftrightarrow \frac{1-z^{-1}}{\Delta T}$$

$1-z^{-1}$ is simply subtracting the past sample from the current sample. To maintain your time axis in seconds, divide by $\Delta T$. If your time axis is in samples (as we often do), the $\Delta T=1$ and further scaling is not needed.

Similarly an integrator which is $1/s$ would be multiplied by $\Delta T$. For all other mappings you see the same $\Delta T$ appear, just as it exists as $dt$ in the continuous time equations.


Finding the impulse response of a system I have the following transfer function.

$$ H(j\omega) = \frac{1+0.5 e^{-j\omega}}{1-1.8 \cos(\frac{\pi}{16}) e^{-j\omega}+0.81 e^{-j2\omega}}$$

I'm trying to find the impulse response of the system. However, I couldn't separate the expression above and I couldn't figure out how I can find the impulse response. Can anybody help me how to solve this equation? Any help would be appreciated.

One can see that the given expression can be decomposed as
$$
\begin{align}
H(j\omega) &= \frac{1+0.5 e^{-j\omega}}{1-1.8 \cos(\frac{\pi}{16}) e^{-j\omega}+0.81 e^{-j2\omega}} \\
&= \frac{A }{1-0.9 e^{j\frac{\pi}{16}} e^{-j\omega}} + \frac{B}{1-0.9 e^{-j\frac{\pi}{16}} e^{-j\omega}}
\end{align}$$

where $A = B^*  = 0.5 - j 3.9375 = 8/16 - j 63/16 $.

Then the impulse respons will be:

$$h[n] = A (0.9 e^{j\frac{\pi}{16}})^n u[n] + A^* (0.9 e^{-j\frac{\pi}{16}})^n u[n]$$

which can be simplified as:

$$h[n] = 2 \cdot 0.9^n |A| \cos( \frac{\pi}{16}n + \angle{A}) u[n]$$

where $\angle{A}$ is the phase angle of $A$. Following is the resulting sequence plotted, from $n=0$ to $n=35$, using MATLAB/OCTAVE.


Calculating an impulse response in a stable environment I would like to estimate the impulse response of a simple environment (room, two speakers, one microphone). I am aware of the issues/limitations of working with "near" LTI systems. That said, I have a rather basic question regarding the mechanics of performing cross-correlation using audio which has been recorded in my environment.

A: my original audio (swept sine) - 1 second in length, 16 bit, 44.1K

B: the result of recording A in my environment

Again, my goal is to calculate the environment's impulse response, h[n], by calculating the cross-correlation of A and B.

Here is my question: I don't know if I should be doing anything to B (trimming on either end, shifting towards the left, etc.) before attempting to estimate the impulse response. The goal will be to use the IR to model what arbitrary audio will "look" like when played in my environment.

Additionally, I'm not sure how long (in # of samples) my IR should be. I've read that using a shorter IR can speed up processing without too much loss in estimating y[n] (the output of the system).

In the final version of this experiment the recording of B will be automated thus I'll need to know how long the recording needs to last to fully capture A's signal.

Thank you in advance.

The easiest way to calculate the impulse response is to divide the spectrum of the measured signal by the spectrum of the excitation and do an inverse FFT of the quotient. That's mathematically equivalent to the cross-correlation method but typically easier. It requires that the excitation signal has a reasonable amount of energy at all frequencies so you don't get any "divide by 0" problems.

The length of the impulse response is mostly determined by the reverberation time of the room. For most of home environments that's well below 0.5 seconds so 1s excitation should be plenty. For convenience I would chose a power of two, which makes the FFTs easier. You can still truncate after the IR has been calculated.

Most acoustic systems are actually strict LTI (if no one moves about), however the signal to noise ratio (SNR) tends to be poor especially at low frequencies. A good way to deal with this is to do "coherent averaging". Create a periodic excitation (just keep repeating the sweep), throw away the first two frame acquisitions (long story) and then keep averaging the acquired frames. 10 averages will give you an SNR gain of 10 dB.

Getting the excitation level right is critical: if it's too low the SNR is bad and if it's too high, something will be overdriven and become non-linear. The speakers are probably the first candidate for distortion.


How can a filter have zero group delay? If you put a wave packet through the passband of a 1st-order low-pass filter, it will be delayed by the group delay of the filter, and remain the same amplitude, right?

If you put the same wave packet through a complementary 1st-order highpass filter with the same cutoff frequency, the group delay curve is the same, so the delay of the packet will be the same, but the gain is much lower, so it will be both delayed and attenuated to negligibility.

Since the output of the highpass filter is very small, if you sum the outputs of these two filters (as in an audio crossover), I would expect it to be negligibly different from the output of the lowpass filter:  Large delayed signal + very small delayed signal = large delayed signal.

Yet if you sum the filter responses, the amplitude is 0 dB everywhere, and the phase is 0 everywhere, and therefore the group delay becomes 0, which would mean that the wave packet comes out with no delay and no changes.  I don't understand how this can be possible.  Don't filters always incur delay?  How can a filter (which also has positive group delay) undo the delay caused by the other channel, especially when this is happening in the stopband?

Which part am I misunderstanding here?

  The best-known crossover types with linear phase are first-order non-inverted crossovers, ... The first-order crossover is minimum phase when its outputs are summed
  normally; it has a flat phase plot at 0°. - The Design of
  Active Crossovers

and

  Here the result of summing the outputs together produces 0° phase
  shift, which is to say that the summed amplitude and phase shift of a
  1st-order crossover is equivalent to a piece of wire. - Linkwitz-Riley Crossovers: A Primer: 1st-Order Crossover Networks

Testing on actual pulses shows how the lowpass (blue) delays the pulse, as expected, and how the highpass (green) can combine with it to produce the original (red) pulse, but how is the highpass pulse occurring before the original if the highpass filter is causal and has positive group delay?  Intuition is failing me.

It does show that the highpass output is not as negligible as I imagined, and the delay is more negligible than I imagined, and as you move the carrier frequency around, these two properties change in a proportional way (smaller delay requires lower amplitude highpass output to correct it).  But I still don't really understand it.

There are a couple of interesting aspects of "reconstruction to unity". First, there are two ways of combining two filters: parallel and in series. For a parallel topology it is ALWAYS possible to find a complimentary filter so that the pairs add to unity. It's easy enough, actually. Simply do $\tilde{H}(\omega) = 1-H(\omega)$. In the time domain that means that the impulse response of the complimentary filter is simply the negative of the original impulse response with 1 added to the first sample. So all the "ringy" stuff cancels out. Now the shape of this complimentary filter is not always what one would expect. For a 1st order low pass it actually is a first order high pass but for higher order filters it tends to have over/under swings in the cutoff region. However, it always exists as a stable causal filter.

Series (or cascade) "reconstruction to unity" is bit more complicated. Obviously the filters would have to be the inverse of each others, i.e. $\tilde{H}(\omega) = \frac{1}{H(\omega)}$. In general this can be done for any minimum phase filter. The inverse of a minimum phase filter is minimum phase as well and both are causal and stable.

So this leaves us with the question of how to interpret group delay in these cases. The cascade case is actually the more interesting. Since the filters are inverse of each other, the phase, and hence the group delay, of one is the negative of the other. So at frequencies were one filter has positive group delay, the other has negative group delay. An easy example would be a low shelf with +6dB of gain and a low shelf with 6dB of cut.
So negative group delays are very real and certainly not a violation of causality. In practice, these show up in areas of the filter that are fairly "non-flat" so the traditional interpretation of "delay of the envelope" doesn't quite apply since there is a fair amount of amplitude distortion as well.

If you Google "negative group delay", you can find a few IEEE articles that have tackled the subject.


How to interpret these MP3 spectrums? Recently, I found I had this song duplicated in my music library, but when listening to them I noticed one version was of noticeably higher quality than the other despite both claiming to have a 320kbps bitrate.

I plotted the spectrograms of some random 5 seconds of each version using python. This is the spectrogram of both channels of the "high quality" MP3:

and this is the "low quality" spectrogram:

And the differences I notice are:

The overall intensity of the low quality version is a lot lower than the high quality file. This difference is more evident as the frequency goes higher.
The low quality song has this weird choppiness above 16 kHz (which I think is the most interesting feature) while the high quality spectrogram seems to be more uniform.

In the power spectrums, I confirmed point number 1 and that the files indeed have a 20 kHz cutoff frequency matching a 320 kbps MP3.

This is the high quality power spectrum:

And this is the low quality power spectrum:

I have been reading about MP3s and first I suspected that the low quality file might be an upscale of a 192kbps (because the spectrum stops being uniform above 16kHz which matches the cutoff frequency of a 192kbps mp3), but I compared my plots to actual 192kbps spectrums and they definitely don't have any information above the 16kHz, so I don't think this is the case.

I also read that this choppiness might be due to the encoding algorithm used, but then that doesn't explain why the power of the whole song is so low (which is easy to notice when you listen to both versions).

Has anyone with more experience with MP3s give me a hint on what may be going on here? I mean, what could cause these differences in both files or any other hints I am missing from the spectrums.

In the end, my goal is to quantify these possible differences to be able to programmatically detect which songs in my library are of "low quality" despite showing a 320 kbps bitrate, so any ideas are also welcomed.

From the power spectrums you have plotted, it seems like there are a lot more deviations with the data. Perhaps when the song was recorded or encoded, the sampling rate (number of samples in each second) was less than the original one.


How to Fourier Transform (basic understanding) I need to fourier transform my light curve (2nd image) into something which looks like the first graph of the first image. I have very little idea about how to actually do this and I only understand the very basics about fourier transforms (the idea, not the mathematics).

T = 1460 #days
z = fft(y)
f = 1/T
plot(f,y)

This is my code attempt so far, but I just get a graph with straight lines or nothing at all. If you could help me with this, it would be fantastic :)

The x axis of my light curve is measured in Julian days and the y axis in flux. When I fourier, I want the x axis to be in micro hertz.

Thank you :)

If you want to plot the magnitude of your data you need to calculate first the absolute value as the output of the fft is complex.

z = fft(y);
z_mag=abs(z); % you could also log-tranform the data for db-scale

To get the right x-Axis you need however the sampling rate. Lets assume that the sampling rate at which the data is recored is f_sampling(some number);

Fs=f_sampling;
n=length(z_mag)
freqHz = (0:1:n-1)*Fs/N;

Keep in mind that this is a double-sided spectrum.


Real-time modulated IIR filter I've been learning to implement an IIR filter for a synthesizer project. It is now more or less stable when the parameters are constant, but ultimately I want the parameters to be arbitrarily modulated.

I figured it would be too slow to compute new recursion coefficients every sample, so instead I'm computing them every "chunk" (a piece of 256 or so samples which the synthesizer produces in a batch and then sends to the audio driver). This brings me to my question:

How should I approach hotswapping the recursion coefficients? Or is there a more correct way of modulating IIR filters?

I am curently linearly interpolating the output of the model generated in the last chunk with the one generated in the current. (at the beginning of the chunk, the output is 100% from the old one, and at the end 100% from the new one) This isn't enough, since the new filter needs to be initialized somehow. I have tried some things with various degrees of success:

The new filter is initialized with signal history pulled from the old one.

This sounds the best, but when too heavy modulation is applied (for example, the cutoff frequency changes too much at once), the new filter might initially output an uncontrollable spike which then peaks. This mostly happens along with a discontinuity in the input signal.

Here's a 60 Hz saw wave LP-filtered using a cutoff oscillating at 20 Hz. (the filter is configured to have two biquads). Such a spike never occurs when the cutoff frequency is kept constant.

The filter history is initialized in a stable state. (e.g. constant signal in and out at the intensity of the last seen sample for LP, or plain zero)

This causes some audible rattle when modulating, no matter how I try to fade in the new filter gradually. The new filter never quite reaches the old one's state.

All ideas and advice are welcome.

Edit:

I have implemented the biquads in a lattice topology, as per Robert's suggestion. Modulation now works great, but raising the cutoff frequency too high (or too low for a high pass filter) breaks something. This should be a single biquad low pass filter with cutoff $\frac{3}{4}\pi$:

The direct form parameters, from which the lattice is derived, work just fine. It is certainly possible, that there is some error in my equations for the conversion. I also noticed that the shape of the frequency response isn't quite as steep as with the direct form. I am suspecting it has something to do with the "tap" coefficients..

What could be going on now? Perhaps I should make a new question with details about my implementation, as the original issue is now more or less solved.

okay, using the snipping tool, i tried to put in just the 2nd-order IIR lattice in the following drawing.  $x[n]$ is being input in the upper left.  the numbering for the $v_*$ coefficients is horrible, but i just cut this drawing out from somewhere else.

now the transfer function i get for this is

$$ H(z) = \frac{Y(z)}{X(z)} = \frac{v_4 + v_3(-k_1 + z^{-1}) + v_2\big(-k_2+k_1(k_2-1)z^{-1} + z^{-2}\big)}{1 + k_1(k_2-1) z^{-1} - k_2 z^{-2}}$$

and if you compare that to the standard direct form

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}} $$

you should have equivalence when

$$\begin{align}
k_2 &= -a_2 \\
k_1 &= -\frac{a_1}{a_2+1} \\
\\
v_2 &= b_2 \\
v_3 &= b_1 - a_1 b_2 \\
v_4 &= b_0 - \frac{a_1}{a_2+1} b_1 + \left( \frac{a_1^2}{a_2+1} - a_2 \right)b_2 \\
\end{align}$$

i would appreciate it if someone else would check my work.  i could well have fucked up a sign or something.

so, wherever you get your five Direct Form $a_k$ and $b_k$ coefficients (e.g., the Cookbook, but remember to divide everything by $a_0$), you can convert them to the five Lattice Form coefficients with the above equations (assuming i didn't mess it up).


Why do linear systems show sinusoidal fidelity? I am looking for a proof for sinusoidal fidelity. In DSP we study a lot about linear systems. Linear systems are homogenous and additive. One more condition it satisifies is that if a signal is a sine or cos wave then the output only changes the phase or amplitude. Why? Why cant the output be a totally different output when a sine wave is given as an input?

A somewhat visual complement to the other answers

You are talking about systems that are linear and time invariant.

Exponential functions have one peculiar property (and can be actually defined by it): doing a time translation results in the same function multiplied by a constant. So

$$ e^{t-t_0}=e^{-t_0}e^t$$

The red exponential could as well be the blue one divided by $e$ or moved 1 second to the right

In general, this also holds for complex exponentials

Can you picture in your mind the plot of an complex harmonic such as $x(t)=e^{j2\pi t}$?
If so, you will see it's like a spring: it rotates along the complex plane as time goes by.

Rotating that spring (multiplying by a complex number in the unit circle) is the same as translating it. You probably have come into this visual effect some time in your life

It's the principle of any standard screw too.

Assume we input this in a linear time-invariant system. You get an output $y$
Now, input a rotated version of this spring. Because of linearity, the output should be $y$ rotated by the same amount.
But since a rotation is equivalent to a time-translation, and the system is time-invariant, the output also has to be $y$ time-translated by the same amount.
So, $y$ has to satisfy the same property as the input: rotating it has to be equivalent to a particular time translation. This only happens when the output is a multiple of the original spring.

How much translation? Well, it's directly proportional to the rotation just like it would happen with a spring. The tighter the loops of the spring (the faster it rotates), the less it time-translates for a certain rotation. The tighter the loops of a screw, the more rounds you have to make it do for it to fit completely. And, when half of the rounds are done, the screw will be half it's way in... The output has to satisfy the same relation so, the output spring $y$ rotates at the same frequency as the input.

At last, a reminder

$$\cos(t)=\frac{e^{j t}+e^{-j t}}{2}$$

$$\sin(t)=\frac{e^{j t}-e^{-j t}}{2 j}$$

So, that thing that happens with exponentials actually doesn't need to happen with cosines and sines in the most general case. But if the system is also real it's a different story...

In general, by this same reasoning, any exponential is an "eigenfunction" (output is proportional to input) of linear time invariant systems. That's why for these systems Z-transforms and Laplace transforms are so useful


Can Principal Component Analysis (PCA) Solve the Cocktail Party Problem? I'm looking into the cocktail party problem and trying to figure out whether something like Principal Component Analysis is enough to separate out all the various voices at the cocktail party into its constituent sound sources.
If its not enough, why? What other techniques should be used in conjunction with it, so that I wind up with distinct signals for each cocktail party patron's voice? Spatial filtering, such as various beamforming methods, has been suggested. But in researching PCA it seems that should (possibly) be enough to "split" the total signal of the cocktail party into individual signals for each human voice attending the party. Beamforming and similiar methods appear to be filters for subsequently focusing on just one of those voices and filtering the rest out.
Can anybody with PCA experience weigh in here with whether or not PCA can solve this, or if its requires additional processing?

The Cocktail Party Problem is a Blind Source Separation (BSS) problem.
Given a linear mixture of signals:
$$ \boldsymbol{y} \left[ n \right] = A \boldsymbol{x} \left[ n \right] $$
We're trying to estimate the signal $ \boldsymbol{x} \left[ n \right] $.
The model can get even more complex with $ A $ being time varying:
$$ \boldsymbol{y} \left[ n \right] = A \left[ n \right] \boldsymbol{x} \left[ n \right] $$
We have 3 main approaches to this problem:

Probabilistic Approach
Looking at the signals as an ensemble of points of a distribution and find the linear coordinate transform to guarantee some property. The PCA approach tries to remove correlation (2nd moment information) while the ICA tries to remove correlation in higher moments (Basically statistical independence).
Time Signal Processing Approach
In case of 2 signals one of them being a reference we can use the adaptive decorrelation filter. Basically we're after removing any time correlation from the signals.
Spatial Signal Processing Approach
We can utilize the known location of the microphones in the room to create adaptive beamforming. The idea is that delayed adaptive summation of the data can change the spatial curve of the array and making matching a certain direction.

Of course in late years we can find work on the subject utilizing Deep Learning approaches. Their main advantage is being able to incorporate additional information (Like using the properties of the signal, be it in a certain language, incorporating visual data on the scene like images and videos [Who is moving the lips when?]).
This is a vast subject and the main idea is to tailor the solution to the specific case of yours.
Modern Robust ICA and IVA (Independent Vector Analysis) can be very effective.
I'd try them first unless you have the case which matches the Adaptive Filter (Which can be proven to be matching the Beam Forming solution under some conditions).


Can you explain about the following formulas? in Signals and System In my textbook, any signals could be expressed like the following formula.

$x(t)=\displaystyle\frac{a_0}{2}+\displaystyle\sum_{n=1}^{\infty}[a_n\cos(2{\pi}nft)+b_n\sin(2{\pi}nft)]$

I understand that any signal could be composed of periodic signals whose frequency is $f$, $2f$, $3f$, $\cdot\cdot\cdot$.

I don't understand the part of "$\displaystyle\frac{a_0}{2}+$".

and the text book said,

$a_n=\displaystyle\frac2T\displaystyle\int_{T_1}^{T_2}x(t)\cos(2{\pi}nft)dt$

$b_n=\displaystyle\frac2T\displaystyle\int_{T_1}^{T_2}x(t)\sin(2{\pi}nft)dt$.

I also don't understand why $a_n$ and $b_n$ are like the formula above.

Thank you for kind answers, and my first language is not English. Please let me read easy English. very Thank you.

Your question is about the Fourier series decomposition of a function (called signal here). First, this does not apply to "any signal" in theory. In practice, it applies to most pratical "periodic" signals. Periodicity is quite important, althought Fourier development may exist for not strictly periodic functions.

In loose terms, it tells you that if a signal is periodic, then inside each period the signal can be expressed as a sum (a linear superposition) of the most simple periodic functions, sines and cosines, with frequencies taken as linear factors of the base frequency.

The $a_0/2$ term comes from the average of your signal over one period $T$. This is due to the fact that sines and cosines do not see constant signals: they are orthogonal to them, since the integral of a sine or a cosine over one period exactly sums to $0$.

The other terms $a_n$ and $b_n$ are related to orthogonality as well, and i guess that you can, with great profit, read the following Wikipedia: Fourier series page in your perferred language.


"Don't window a transient signal" I'm a SP student and still a beginner in the field.

When we were learning about window functions, my professor said "you're gonna almost always want to window your data but don't window a transient signal." It's a (semi-)online course with lecture videos so I couldn't ask for an explanation.

Can anyone explain what a transient is and why we wouldn't window it?

To build on Laurent's answer, here is an example.  Top frame shows an example transient signal: a damped sine wave.  As the signal decays very quickly, the first 0.1s of the signal is the most interesting.  Second frame shows hann window.  Hann window is almost zero over that first 0.1s.  Third frame shows what happens when you apply hann window to the signal, the most interesting part of it is completely thrown away.


If an autocorrelation function of a process is a rectangular function, then is the process deterministic? I would like to clarify one doubt regarding autocorrelation function $R_x (\tau)$ of a process $X(t)$. Suppose $R_x (\tau) = \Pi (\tau)$, where $\Pi (\tau)$ is the rectangular function. Then is the process $X(t)$ deterministic.

An autocorrelation function like that is not possible for ordinary signals. The autocorrelation function is the Fourier transform of the power spectral density, which is a strictly positive quantity. However, to get a rectangular function as Fourier transform you need to have negative values. That is in conflict with the strictly positive power spectrum.


What is this sinc doing in my LP-turned-HP filter output? I'm trying to HP-filter a signal by LP-filtering it and subtracting the output from the original filter. To obtain the LP-filtered signal, I backward-forward filtered using a third-order Butterworth filter for a 16,000-sample signal. The cutoff frequency was 0.5 Hz @ 250 Hz. No window was used.

In the first 85% or so of samples, the output of this method seems to overlap completely with the actually HP-filtered signal. However, at the end of the output, there is a supsiciously sinc-like phenomenon.

What causes it?
What could I change to make this method feasible?

Here are the two signals overlaid. The blue one is HP-filtered, and the red one obtained by subtraction.

Here is the absolute difference per sample. They are never exactly equal; the minimum difference is 8e-12. 87% of the values are below 0.01.

(If there's a better way to obtain an HP filter from an LP filter's $H(z)$, I'd be happy to learn it.)

You are doing a backwards-forwards filter for one, and a forwards-only filter for the other.  Thus the two filters will behave differently at one or both edges.

If your signal is strongly stationary and exactly periodic in the data window length, you could try wrapping the data around circularly and start your filters outset by their impulse response lengths.


Signal Processing and Emotion Recognition I have to write a paper on emotion recognition.
More specifically on how different parameters affect the algorithms effectiveness. I have already read a lot about this and done research using OPENSmlie, (extracting mostly MFCCs) and WEKAs clasificators to get the best results. I use EMO DB.

However what i need is more theoretical knowledge, and information about already published papers on that subject. Especially newer ones, say from the last couple of years. I want to know, how peaople approach this problem, and how good their latest results are.

I was wondering if maybe there were some experts on that here.
Could you recommend any good recources / links to such papers. Or an effective way of finding the latest ones?

For the latest improvements read the papers from the latest speech conferences, for example here are two tracks on Emotions in latest ICASSP:

Emotion recognition I

Emotion recognition II

You can also check reviews:

Speech emotion recognition: Features and classification models

Emotion recognition from speech: a review

Survey on speech emotion recognition: Features, classification schemes, and databases

To find the relevant research use Google Scholar.


GPS Receiver in Frequency Domain I'm trying to transmit and receiver the GPS payload in both Time and Frequency Domain.

For transmitting and receiving in time domain :

numBits = 10;                        % LFSR length
codeLength = 2^numBits-1;           % Spreading code length

payload1 = 1*(rand(1,100)>0.5);
payload2 = 1*(rand(1,100)>0.5);
payload3 = 1*(rand(1,100)>0.5);

payload1 = 2*payload1 - 1;            %BPSK Modulation
payload2 = 2*payload2 - 1;            %BPSK Modulation
payload3 = 2*payload3 - 1;            %BPSK Modulation

%  payload1 = ones(1,10);
%  payload2 = ones(1,10)*2;
%  payload3 = ones(1,10)*3;
%% PN code generation

% initial sequence is taken from IRNSS document
goldCode1 = generateGoldCode(935,239,codeLength,numBits);
goldCode1 = 2*goldCode1-1;                        % Assign code to 1 to 1 and -1 to 0

goldCode2 = generateGoldCode(38,381,codeLength,numBits);
goldCode2 = 2*goldCode2-1;

goldCode3 = generateGoldCode(564,561,codeLength,numBits);
goldCode3 = 2*goldCode3-1;

%% At transmitter part
payloadData1 = kron(payload1,goldCode1);         % Product of payload data and gold code to spread the sequence
payloadData2 = kron(payload2,goldCode2);
payloadData3 = kron(payload3,goldCode3);

TxData = (payloadData1 + payloadData2 + payloadData3);

%% At receiver part
%% Time Domain
crossCorr = xcorr2(TxData,goldCode2);
outputData = crossCorr(codeLength:codeLength:end)/codeLength;

For Frequency domain I want to receive each bit at a time:

starting =1;
i=0;
crossCorr_temp=[zeros(1,100)] ;
goldCodeFFT = conj(fft(goldCode2));
for val = 1023:1023:length(TxData)
  payloadFFT = fft(TxData(starting:val));
  crossCorr_FFT = ifft(payloadFFT .*goldCodeFFT);
  corr_time = crossCorr_FFT(1:end)/codeLength; **#Doubt here**
  starting = starting +1023;
  i=i+1;
end

How do I approach this frequency-domain way to receive the GPS payload bits?

I do not see any issues with your approach for the frequency domain computation, given that the following is a circular convolution in the time domain of the two signals a and b:

$$corr(a,b) = ifft(fft(a)fft(b)^*)$$

But to receive each bit you need to evaluate the maximum in each interval to determine if it is positive or negative.

Further in actual application, the received signal would be sampled higher than one sample per bit, with the unknown Doppler offset frequency removed (carrier recovery) and aligned such that your start of the FFT is at the start of the symbol. There are also 20 1023 length symbols per actual data bit in GPS.

In this case the FFT would be rate matched to the received signal. You could also process longer blocks as I show in the example plot below, by simply zero padding the rate-matched prn sequence out to the length of the data block prior to taking the FFT.  In the plot below this was done for 1e6 samples of actual received signals sampled at 20 MSps.


Fourier transform of an image I am doing a project where I need to take the FFT and IFFT of a photo in MATLAB. There is a principle in optics where the intensity distribution from a lens is equal to the Fourier transform of the aperture that the light initially goes through. I'm having issues finding the FFT and IFFT of an ideal source.

A circular aperture should have an intensity pattern called Airy disk, which should result from the Fourier Transform of a circle. Here is a photo of the Airy disk that I'm using in my code:

Taking the inverse Fourier transform of the Airy disk should result in an image of a circular aperture, but all I'm seeing is black when I convert to uint8. If I leave it as type double then I see things that also look like Airy disks in the corners. Here is the code I'm using:

pic=imread('airy.jpg');
pic=rgb2gray(pic); %Convert to grayscale
imshow(pic)
title('Initial');

pic=double(pic); %Convert to double before transforms
pic=ifft2(pic); %Take IFFT of Airy disk
pic=real(pic); %Take real part

figure()
pic=uint8(pic); %If commented out, black screen has Airy disks in corners
imshow(pic);

Does anyone have any ideas on why the IFFT of the Airy disk doesn't return a circle? Thank you

You have the gist of it, but there are a couple of problems. I'm sorry, I don't have MATLAB, but I'll do my best to give you MATLAB-like pseudo-code to work with.

The main problem is that your Airy disk is an image, i.e. amplitude only (the result of taking abs). To reconstruct the circular aperture, you need phase as well. Here's code to first generate your complex Airy disk (PSF), from which you can then reconstruct the circular aperture (pupil function). Starting with this file, circ.jpg:

pupil = imread('circ.jpg');
pupil = rgb2gray(pupil);
% use fftshift to block-swap the fft output (move DC to the center)
complex_airy = fftshift(fft2(pupil));

% now, starting with this complex airy pattern, you were close:
pupil_reconstructed = ifft2(complex_airy);
pupil_reconstructed = abs(pupil_reconstructed);

Below are the amplitude and phase of the complex_airy pattern, obtained using abs(complex_airy) and angle(complex_airy), respectively. (I've included only the centers of the images, to show detail.) The problem with your initial approach is that it's an ill-formed inverse problem: a single Airy amplitude pattern can correspond to an infinite number of pupil functions. When you inverse transformed the amplitude only, you assume the phase is everywhere zero, which--as is clear from the phase image below--it's not. There are statistical approaches for estimating the pupil function from the Airy amplitude alone, using lots of other constraints gleaned from the light source properties and optical system, but I don't think that's what you were after.


How can I reduce background noise from recorded audio conversation and separate speakers in python? I have a recording of an audio conversation of two persons. I want to remove the background noise from audio and then also separate both persons' audio separately in python.

I have gone through some concepts about applying filters on the audio, but I was unable to implement the correct filter. In some audio recordings, it able to cancel out the noise partially.

For Separation of speakers, I found some libraries in python like PyAudioAnalysis, uis-rnn (Speaker-diarization), but they have their own trained models, which are solving the problem. Any paid APIs are available to solve this speaker-diarization problem?

Is there any way to solve this problem?

Thank you:)

I will contribute to this question from my area of experience, which is in speech enhancement. For noise suppression, you can use classical speech enhancement algorithms (spectral subtraction, wiener filtering, etc) but they mostly introduce artefacts in the enhanced speech. Additionally, there are also machine learning approaches which you can try. I am providing you with one of the python-based denoising libraries here:
https://github.com/staplesinLA/denoising_DIHARD18

Remember, for commercial usage you need to ask the author and the group.


Noise Removal from an Image Using OpenCV (Comparison to Neat Image) I tried removing noise from the image shown below using Median Blur in OpenCV. But i'm not able to remove the colour noise completely as it is done in Neat Image. Any suggestions.?
 1. Original Input Image

 Median Blur Output
 Neat Image Output

NeatImage probably uses Wavelets based Noise Reduction.
You can look for methods based on that.

Today you need methods which are "Edge Aware", namely they smooth yet keep edges in tact.

Have a look here:
https://github.com/RoyiAvital/Fast-Anisotropic-Curvature-Preserving-Smoothing


The way to measure chromatic tuners' precision I have a question about measuring precision of chromatic tuners.

I want to divide 2 cases.

CASE1:

link; Korg AW-2G's precision

If you check the link or this picture, it says that
Korg AW-2G has the precision of +/- 1 cent.

What's the way or measurement to get this result (+/- 1 cent)?

CASE2:

There are people using Arduino (ATMEL328 MCU) with microphone or piezoelectric diaphragm for the ADC input.

One blog introduced the autocorrelation method
for detecting the instrument's pitch.
After calculating autocorrelation, using parabolic interpolation gives
"accurate" results. (This is not sarcasm, don't get me wrong. :D)

Suppose I made a chromatic tuner using ATMEL328 and a piezoelectric diaphragm.
Also, suppose I chose sampling rate 9240 Hz, 512 samples, the autocorrelation method (with parabolic interpolation).

In this case, how should I define the precision of this chromatic tuner?
I can't buy tuning forks to compare the results.
Let's suppose I'm testing with electrical guitar, classical guitar, and violin.
So any ideas?

Best Regards~

Precision and accuracy are 2 different things.

1 cent is 1/100th of the frequency difference between 2 adjacent semitone pitches in equal temperament tuning.  The frequency ratio is exp(ln(2)/1200) or roughly 1.00059

Precision of 1 cent might mean the tuner reports frequency error from the nearest semitone in 1 cent units, whether accurate or not.

To test accuracy, you might want to synthesize test signals that "sound like" your desired instrument (in software, or using a calibrated programmable synthesizer).  Feed those test signal to your tuner (or your tuner's pitch estimation algorithm in a software environment as unit testing).  If the results matches the synthesized pitch to within 0.5 cents, then your accuracy is 1 cent at that frequency.  Vary the test signal (perhaps randomized) to different semitone notes and cent (or cent fraction) offsets over your tuners entire range to find the worst case to evaluate (get a statistical estimate of) your tuners accuracy.

You can add Gaussian random white noise or recorded typical background room noises to your generated test signals to test measurement precision (repeatability) at realistic S/N ratios as well.

One "gold standard" is to measure the WWV audio frequencies (500 and 600 Hz), and use those to calibrate your synth and test your tuner.


Positive slope of BER with a greater frequency deviation? I am trying to implement a frequency-modulated communication channel on MatLab. However, there is something I don't understand about the sensitivity of the frequency modulator (or the frequency deviation as my signal has an amplitude equal to 1). My message has a 500Hz frequency and is sampled at 20kHz. I added white gaussian noise to my bandpass signal (carrier of 6kHz).

Why does the BER of the demodulated message have a positive slope after 250Hz ?

NB : it is written [dB] for the BER on the graph but it is not, I forgot to remove it

This could happen as discriminator gain is increased with a filter discriminator approach since in many of those approaches the gain would be maximum and linear for small signals only and then the slope of the discriminator slowly goes down coinciding with the results in your plot (such that you no longer get a perfect sine wave out for a sine wave in—- so you can also look at it as converting some of your signal components to other harmonics at the discriminator output due to that non-linearity). Ultimately as you increase the gain further your signal will start to saturate and then even invert depending on the wider bandwidth shape of your discriminator.

Plot the derivative of your frequency discriminator to see this more clearly as to its usable range for a given frequency deviation- the ideal for no distortion is to have a constant slope.

I see from the link you provided that the slope is indeed constant over a limited usable range. So specifically in your case I assume as $k_f$ was increased, it would cause a proportionate amount of the signal to go beyond the discriminator range (where your response is 0).


What does mean and variance do in Gaussian noise? I = imread('us-currency.png');
mean = 0.5;
variance = 0.5;
J = imnoise(I,'gaussian', mean, variance);
imshow(J);

What are the effects of mean and variance on that image?

Does mean increases the brightness?
Does variance increases the amount of noise?

Input

Output

Does mean increases the brightness?

As you add a constant mean to each pixels of the image you will get an image that has an overall higher average image pixel value. Since images are generally encoded with the maximum pixel value being white, you would get an image which is closer to white, and hence perceptually appears brighter.

  Does variance increases the amount of noise?

That one deserve a little bit more background. Indeed, quoting from Wikipedia:

  In experimental sciences, noise can refer to any random fluctuations of data that hinders perception of an expected signal.

So in that context, even adding or subtracting a constant value to the entire image could be coined as noise (especially if this makes the image too bright or too dark so as to make it harder to discern objects in the image). If on the other hand we focus our attention to what is commonly described as a "noise" effect then this typically refer to some kind of perceived "dots" that occurs due to somewhat large variations in pixel values in close proximity. Since the variance is directly linked to the variation in pixel values, it is clear that an increase of the variance would increase the amount of noise (in the sense of perceived quantity of "dots").


Difference equation system identification Suppose I have a difference equation modeling a vehicle like this:
$$d[k+1]=d[k]+a\cdot u[k]+b,\tag1\label{eq}$$
where $d[k]$ is total distance traveled at time $k$, $u[k]$ is engine input at time $k$ (e.g. some measure of engine exertion at that time, not really important exactly what it is), and $a, b$ are parameters that I want to estimate from data points that I measure, for example by least squares.
From the difference equation I see that the units of $a\cdot u[k]$ and $b$ need to be distance. But suppose that the data I gather are actually samples of the vehicle's velocity (not position) for some range of input values, so I can plot velocity vs. engine input.
My question: I'm confused because it seems that it wouldn't make sense to do least squares regression to fit $a\cdot u[k]+b$ to these points since that would violate the units, but that is what is done in an example in my course. After the parameters are determined, the example then uses those same parameters to model the distance traveled over time via $\eqref{eq}$, going off of some initial starting distance $d[0]$. Am I missing something in this example?

Note that for a small sampling interval $T$, $\big(d[k+1]-d[k]\big)/T$ is a good approximation for the velocity. So if you fit $au[k]+b$ to a given set of measurements $v[k]$, it is valid to conclude
$$d[k+1]=d[k]+T\big(au[k]+b\big)\tag{1}$$
In the text you refer to they might have normalized $T$, so it changes the units without changing the values of $a$ and $b$.