texts: Keysight ADS Gaussian noise + understand signal plan (BPSK-Mapping a. Demapping)

texts: I have already created the circuit in the Schematic.I am concerned with the understanding of the diagram. In the diagram you can see the Gaussian noise before and after the addition and the question is why before the addition and after the addition both signals look the same and how you have to understand the output.


texts: Which transform most closely mimics the human auditory system?

texts:
The Fourier transform is commonly used for frequency analysis of sounds. However, it has some disadvantages when it comes to analyzing the human perception of sound. For example, its frequency bins are linear, whereas the human ear responds to frequency logarithmically, not linearly.

Wavelet transforms can modify the resolution for different frequency ranges, unlike the Fourier transform. The wavelet transform’s properties allow large temporal supports for lower frequencies while maintaining short temporal widths for higher frequencies.

The Morlet wavelet is closely related to human perception of hearing. It can be applied to music transcription and produces very accurate results that are not possible using Fourier transform techniques. It is capable of capturing short bursts of repeating and alternating music notes with a clear start and end time for each note.

The constant-Q transform (closely related to the Morlet wavelet transform) is also well suited to musical data. As the output of the transform is effectively amplitude/phase against log frequency, fewer spectral bins are required to cover a given range effectively, and this proves useful when frequencies span several octaves.

The transform exhibits a reduction in frequency resolution with higher frequency bins, which is desirable for auditory applications. It mirrors the human auditory system, whereby at lower-frequencies spectral resolution is better, whereas temporal resolution improves at higher frequencies.

My question is this: Are there other transforms which closely mimic the human auditory system? Has anyone attempted to design a transform that anatomically/neurologically matches the human auditory system as closely as possible?

For example, it is known that human ears have a logarithmic response to sound intensity. It is also known that equal-loudness contours vary not only with intensity, but with the spacing in frequency of spectral components. Sounds containing spectral components in many critical bands are perceived as louder even if the total sound pressure remains constant.

Finally, the human ear has a frequency-dependent limited temporal resolution. Perhaps this could be taken into account as well.


texts: Z domain transfer function to difference equation

texts:
I want to convert this transfer function: $$\ \frac{2\cdot(z-0.5)\cdot(z-0.6)}{z-1}$$ to a difference equation, but I have no idea how.

Can someone help me?

Thanks,

Arjon


texts:
What is the difference between soft thresholding and hard thresholding. Where we use soft and hard thresholding in image for denoising. I understand that in hard thresholding, the coefficients below threshold value are set to zero and the value above the threshold is set to one. Please explain me about soft threshold.
Please explain whether the threshold value is the intensity value of the image. For example if the intensity value ranges between 0 to 255. In case of hard thresholding if the threshold value is considered as 100 then the values below 100 is set to 0.The value above 100 are retained. Is this correct? Please correct me.


texts: Using MATLAB to plot the input and the magnitude spectrum of the signal

texts:
I have an aperiodic signal $v_{out} = e^{-t} u(t)$ (real exponential signal) from discharging capacitor.

I was trying to plot using MATLAB 15 seconds of this signal in time domain? I am thinking how to plot because I donot know $u(t)$. But I am seeing many sources they involve the term $u(t)$ and the Fourier transform is independent of $u(t)$.

I found the Fourier transform  $V_{out}$ of $v_{out}$ to be $\frac{1}{1+j\omega}$ where $\omega = \frac{2\pi}{T}$.

$|V_{out}(j\omega)| = \frac{1}{1 + \omega^2} $

How do I plot the magnitude spectrum?(Is it the plot betwen $|V_{out}(j\omega)|$ and $\omega$?)


texts:
I need to implement an encoder which compresses a 5-dimensional structure of 10 bits values. Each dimension has between 4 and 12 elements. If a dimension ever has more than 12 elements, it is partitioned in half.
So far, I have been using a separable DCT to do that, but my implementation is very slow: I have implemented a naive DCT transform, which has complexity of $O(N^2)$.
I have tried some fast algorithms, but they seem to be slower than the naive approach:
a naive DCT for 6 elements is faster the Fast version, but for $N = 4$ and $8$, the fast approach is definitely faster. Also, I cannot force the friendly power-of-2 sizes.
I remember the Cooley–Tukey FFT algorithm decomposes a sequence of size $N$ into their factors $N_1$ and $N_2$ such that $N = N_1N_2$. I also know that there are specialized algorithms for 2D DCT due to their use in video and image codecs.
My question is: Is there any reference about a fast multidimensional DCT or even a "reversed" Cooley-Tukey
for DCT of which I could transform my 5D array into 1D and use a fast DCT?
Exrta bits:  We are researching not only the DCT-II but the whole family of trigonometric transforms., insights about any of the DST and DCT are welcome!
Edited 9/23
Just as an extra comment, this solution is being implemented in C++.


texts: Derivative of signal with missing samples

texts:
I have software that tracks an object moving (in the x-dimension only) across a video shot from a stationary camera. I need to find the velocity and acceleration of the object as functions of time. This calculation can be done offline, so there are no real-time constraints. The motion of the object is also pretty regular, decelerating approximately constantly throughout the duration of motion.

Normally I would design an FIR differentiator to do the job, however the object in the video becomes occluded for up to a few frames at a time on a regular basis throughout it's motion, so I don't have a point for it for several frames. In other words, my data looks like this:

frame_index = [1 2 3 6 7 8 9 13 14 15 16 17 18 21 22 ... ]
x_position  = [0 0.2 0.39 1.14 1.22 ... ]

I've thought about interpolating and then using an FIR differentiator but am unsure what a good scheme for interpolating irregular data like this would be. I've also thought about using a formulation based on Lagrange interpolating polynomials to calculate derivatives directly, but do not know enough about them to understand trade-offs. My goal is accuracy.

So, what would be the best way to go about calculating approximations of the first and second derivatives of the data?


texts: How to come true the convolution with a 2-D array in frequency domain with python

texts:
I use the pretty simple example used in many books to understand the convolution in the frequency domain. But I cannot find the real results.
import numpy as np

img = np.array([[0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 1, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0]])
kernel = np.array([[0, 0, 0, 0, 0],
                   [0, 1, 2, 3, 0],
                   [0, 4, 5, 6, 0],
                   [0, 7, 8, 9, 0],
                   [0, 0, 0, 0, 0]])

fft_img = np.fft.fft2(img)
fft_kernel = np.fft.fft2(kernel)

fshift = (fft_img* fft_kernel)  # Convelution theory at frequency domain

img_back = np.fft.ifft2(fshift)

img_back = np.real(img_back)
img_back = img_back.astype(np.int)

print(img_back)

The output of the img_back is
[[9 0 0 7 8]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [3 0 0 0 1]
 [6 0 0 4 5]]

However, the answer should be
[[0 0 0 0 0]
 [0 9 8 7 0]
 [0 6 5 4 0]
 [0 3 2 1 0]
 [0 0 0 0 0]]

So, is there someone who can help me? I don't know the reason why. And how to come true the ideal result?


texts: Is there a way to extract the difference between the 2 signals? Fitting?

texts:
I'm very new to signal processing and have been struggling this far.
I have a device that measures the signal due to the magnetization of iron particles.
In order to calculate the response of the particles, I did 2 measurements:

Background measurement, so without particles, so the signal should be 0, but as the figure shows, it is not.
Measurement with particles.

The problem is that even at the background measurement I get a sort of signal, because of many complex reasons like field through, in the device, there is iron because of the production process and so on...
It is true that the magnitude of the signal of the base measurement is 2 orders lower than the one with the particles but the shape of the signal and the peaks are already to be seen at the base measurement. So one should say, that the signal I got while having the particles in the device, does not only come from the particles.

Question: now my question is how to extract the signal that only comes from the particles? A simple subtraction is not the answer, as the magnitude is not alike. Filtering in the frequency domain did not help as the particles are simulated with a frequency f0 and in the signal, we just collect the harmonics.
Can please anyone help?

My suggestion:

take out the peaks at 44 and 65 because they are certainly not from the particles, as one can definitely see them in the base measurement without the particles.


texts: How to find the impulse response of filter that outputs the autocorrelation of the input?

texts:
I am trying to solve question 9.61 in the attached image.

I guess, by evaluating the autocorrelation for a specific instance,  say $ t = 0 $ , that the impulse response, $h(t)$, should be : $$ h (t) =  x(-t) $$.

But can I get this from equating the convolution integral to the given output integral, i.e.

$$y(t) = \int_{-\infty}^{\infty} x(\tau)x(t + \tau) d\tau = \int_{-\infty}^{\infty} x(\tau)h(t - \tau) d\tau $$ .

Is it then right to equate as below?

$$ x(t + \tau) = h (t - \tau) $$

Thank you!


texts:
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.


texts:
I have to find the range of r which makes H(z) stable.  There is no restriction of left sided or right sided.
Then, both r<1 (z>r) and r > 1  can be the answer?


texts: matlab code question with fractional fourier tranform

texts:
I want to estimate the parameters of a chirp signal using the fractional fourier transform. Following the paper estimating chirp parameters, the chirp rate $\mu$ can be obtained with (minus) the cotangent of the rotated angle $\theta$, $\mu=-\cot{\theta}$. The simulation code I wrote is
   fs=256;
   t = 0:1/fs:1023/fs;
   y = exp(-j*24*pi*t.^2);
   XF = [];
   for alpha = 0:0.001:1
        XF = [XF,frft(y,alpha)] ;
   end;

   figure;
   surf((abs(XF)).^2);
   colormap('jet'); view(0,90);
   shading interp;

The discrete fractional fourier transform function frft is adopted from code.
The figure peaks at $\alpha = 0.772$. As a result, the rotated angle is $\theta=0.772*\pi/2$ and the estimated chirp rate is $\cot\theta=0.3743$. Obviously, this is not the desired result.
What may be the problem with my code?  Thank you!


texts: Finding the parameters and endpoints of a sinusoid

texts:
Suppose I have some signal

$$
s(t) = n(t) + \left\{\begin{aligned}
&0 &&: t < t_0\\
&A e^{i (2 \pi f t + \theta_0)} &&: t_0 \le t \le t_1\\
&0 &&: t > t_1
\end{aligned}
\right.$$

with unknown beginning and end times, $t_0$, $t_1$, and unknown frequency, initial phase and amplitude $f$, $\theta_0$, $A$. The signal also has some complex Gaussian white noise $n$ with a variance of $\sigma_n^2$, giving this signal an SNR of $A^2/\sigma_n^2$.

My typical SNR is about 0.0025

I am interested in finding the unknown parameters above, particularly the beginning and end times to as high precision as possible.

I am also interested in knowing the theoretical limit on the precision of the start and end times vs SNR and signal length/number of samples. As the SNR increases information which determines $f$, $t_0$, and $t_1$ is lost. It seems like there should be some abstract, information theoretic way to relate the SNR to the estimation accuracy of these variables which is method agnostic.

If this is possible I can figure out whether or not I can reach the required precision on these variables by improving my analysis software or if I need to make hardware changes to improve the SNR.

My current method is to use an FFT to identify $f$, and then bandpass around it to remove most of the noise. I compute the amplitude of the signal. The noise amplitude distribution is second order chi ($\chi$) where $\sigma_n$ is the scale of the distribution. I fit the noise in my band passed region to this distribution by maximizing a likelihood function on $\sigma_n$ over a region where I know there is no signal.

$$\mathcal{L}(\sigma_n) = \sum_\text{noise samples} p_\chi(\sigma_n)$$

Then I maximize a likelihood as a function of $t_0$, and $t_1$ like so:

$$\mathcal{L}(t_0, t_1) = \sum_{t=0}^{t_0} p_\chi(\sigma_n) + \sum_{t=t_0}^{t_1} p_\text{norm}(\mu_s, \sigma_s) + \sum_{t=t_1}^{t_\text{max}} p_\chi(\sigma_n)$$

The signal distribution is a non central chi distribution, but that approximates to a Gaussian distribution reasonably well in the range I am working in. I am setting $\mu_s$ and $\sigma_s$ to the sample mean and standard deviation for the $t$s between $t_0$ and $t_1$ so they don't need to be optimized over.

This method gives me $f$, $t_0$, and $t_1$ reasonably quickly accurately, but is there a way to do it better?


texts:
I am a beginer in signal processing.
I am currently studying a signal $s(t) = \cos(h(t))$, and I have to calculate its autocorellation.
Since $h(t)$ is periodic of period $T_s$, I calculated the value of

$$R_s(\tau) = \frac{1}{T_s} \int_0^{T_s} \cos(h(t)) \cos(h(t-\tau)) dt$$

Honestly, this calculation is not easy (since $h(t)$ is also a cosinus), and in order to simplify my final document, I would like to use complex representation by writing $s(t) = Re\left\{ e^{ih(t)} \right\}$. However, This cannot work for the autocorrelation since $$Re\left\{ e^{ih(t)} e^{-ih(t-\tau)}  \right\} = s(t)s(t-\tau) + Re\left\{ \sin(h(t)) \sin(h(t-\tau))  \right\}$$

So, do you know a way to calculate $R_s$ by calculating $\frac{1}{T_s} \int_0^{T_s} e^{ih(t)} e^{-ih(t-\tau)} dt$ ?

I was thinking about a relation like $R_s(\tau) = \frac{1}{2} Re \left\{\frac{1}{T_s} \int_0^{T_s} e^{ih(t)} e^{-ih(t-\tau)} dt \right\}$, but I cannot prove it :(


texts:
signal processing is not my domain of expertise , but I need to make use of it, to do an HRV analysis. These are steps I do starting from my raw ECG signal.

ECG signal -> RR peak detection -> find RR intervals -> cubic spline interpolation -> IIR high pass filter -> FFT

The peaks I get is a mirrored image, I guess its because its a pure signal from a digital patient simulator like all at 60 bpm or all at 120 bpm etc and follows the rule of Hermetian symmetry.

However my FFT output at pure 60 bpm should have indicated a peak at 1 Hz and a 120bpm at 2 Hz, is what I understood. But my all peaks start from 0.

Please help

PS : Implementing this all in an android application.

PS : This is what my output looks on a 60 bpm pure signal.

PS : Added log data

FFT output : http://www.yourfilelink.com/get.php?fid=1337132


texts: Balance brightness and contrast between two images

texts:
I have two different images for the same object. Each one was taken from different device (scanner, mobile camera, professional camera ..etc).
the process of  registeration one of them to the other was done correctly. Now I need to do a subtraction between to them for comparison purpose. The problem that is the contrast and brightness are pretty different. I tried these two methods (I tried them on the whole image(globally) and on pairs of the two images (locally)):

Gray-Scale method:

ImageB= ImageA * (Mean_Of(ImageA)/Mean_Of(ImageB));
Color method:
Convert it to HSV then:

ImageB[S]= ImageA[S] * (Mean_Of(ImageA[S])/Mean_Of(ImageB[S]));
ImageB[V]= ImageA[V] * (Mean_Of(ImageA[V])/Mean_Of(ImageB[V]));

return to RGB again

However, the two methods did not give me the results I want. They are good and work well on many images but not in all cases.

Any suggestions?
Many thanks


texts: Camera calibration and extrinsic parameters for perspective transformation

texts:
I'm working in a photogrammetry camera and laser application where I need to measure a displacement in a specific plane, this is a quick diagram of the setup:

In order to get proper measurements I'm trying to do a projection from the 3d points to the 2d camera plane  proceeding like this:

First I'm doing the Camera calibration using opencv and a chessboard,
I'm taking a few chessboard shots in different angles and applying the
function initUndistortRectifyMap, I´ll have the distortion coefficients,
intrinsic and extrinsic parameters.
With the intrinsic parameters and the coefficients I'll apply lens correction.
The extrinsic parameters of the previous step are for each chessboard
position and I assume that they won't be valid for my particular
application, so I'm placing a chessboard in the above marked area
as "Measurement Plane" and running again the initUndistortRectifyMap function
which this time it's giving me, I think, correct extrinsic parameters
of the ROI.
Once I've got the "correct" intrinsic and extrinsic parameters I´m
applying the camera matrix to calculate the 2d points (u,v) projection:

My questions are:

Is it the correct way to proceed for this application?

The camera system only reads X,Y pixel points each time, but to perform a 2d projection seems like I'll need the Z, how could I have this last coordinate of the 3d point?


texts: Fourier Transform minimum sample requirement

texts: Assuming a signal is sampled adequately, what is the minimum size of an FFT window that would allow detecting a specific frequency?  Is it necessary to have samples for at least one complete period at that frequency?


texts:
I have to make comparison between 155 image feature vectors. Every feature vector has got  5 features.
My image are divided in 10 classes.
Unfortunately for using support vector machine i need at least 100 images for class, There is any alternative?


texts: Stochastic process inference from partial observations

texts:
Consider a set $U$. My signal is a piece-wise constant "function"
$Sig: t \mapsto s$, i.e. the signal at time $t$ equals to some subset
$s \subset U$. One can see $Sig(t)$ as a stochastic process.

For a given sequence of time points $\{t_i\}$, we have sequence of
partial observations $\{ Ob(t_i) \}$. Each observation $Ob(t_i)$ is
not a noised value of $Sig(t_t)$, but rather a "random" part of the "real"
signal. I mean:

$$ Ob(t_i) \subset Sig(t_i) \,.$$

(For example, $Sig(t_i)$ can be a graph, and $Ob(t_i)$ will result in a random subgraph of $Sig(t_i)$).

Is there a theory that describes such partially observed stochastic
processes?  In particular, is there a method to infer (estimate) the
parameters of the underlying stochastic process $Sig(t)$ using only a
given sequence of observations $\{Ob(t_i)\}$?

If you want to add some addition properties (constraints) to $Sig(t)$ or to $Ob(t)$, you are welcome.

Update:

For example, a very simple problem from this class is considered here: https://stats.stackexchange.com/questions/83998/if-maria-performs-more-observations-per-unit-of-time-than-maximilien-how-can-he


texts: When do phases not exist for spectrograms?

texts:
I have been reading a paper on the "Single pass spectrogram inversion"

and I came across this in the Introduction part.

  In many applications, the analysis and modification of the Short-Time
  Fourier Transform (STFT) and the Short-Time Fourier Transform
  Magnitude (STFTM) of speech and audiosignals are necessary. These
  applications include, but are not limited to audio enhancement,
  reverberation analysis, time and pitch modification, and noise
  cancellation. Phases are either lost, become meaningless as the
  spectral representations are manipulated, or simply do not exist for
  artificially constructed spectrograms. The objective, then, is to use
  these spectral representations to generate a real-valued signal that
  corresponds as closely as possible to the original spectrograms.

I do not understand why phases do not exist for artificially constructed spectrograms. A spectrogram simply represents what frequencies exist for a certain duration of time, right? So technically, we could just have zero phase for all the frequencies involved and it would be a realizable signal.

Any help is appreciated


texts: How to properly detect segments in a LCD segment display

texts:
I have been given the task to develop a code which checks if the display of a device is showing information correctly (i.e. any segments are activating when they are not supposed to).

To meet this goal the device is able to show different segment sets, especifically symbol segments, numeric segments, all segments and none.

The image is captured in these conditions:

To prepare the image, I rotate it and  afterwards apply some filtering while decomposing the image channels and working on the Green channel, resulting in the following image.

I am using HALCON to carry out this task, developing a code which performs a pattern matching over a set of ROIs with the following result. For the following images I zommed in to the region most interesting for the pattern matching.

As you can see in the resulting pattern, there are some segments I am not checking (i.e the decimal points) while some not being completely checked like the PROG segment.

After altering the parameters I get the following pattern:

As you can see, it now checks for more information, but now there have appeared undesired traces in the pattern, which might end up affecting the performance of the pattern matching procedure.

Note that the images taken are not the original ones, I have performed the following filters on the image:

shock_filter (ImageRotate, SharpenedImage, 0.1, 20, 'canny', 1)

binomial_filter (SharpenedImage, ImageBinomial, 9, 17)

My questions are:

Am I approaching this pattern detection correctly? If not, how should I approach it.
For better results, should I illuminate the display as to reduce the amount of noise to a minimum?

I feel like I should do a matching per every segment group (Perform a matching for every ROI)


texts: How do I examine if the signal can be z-transformed?

texts:
It's given signal

$x[n]=\sin(\frac{2 \pi}{N} m n) u[n]$

where $u[n]$ is the unit step function. Can I calculate Z-Transform?

$\mathcal{Z}$ transform exists when
$$ \sum_{n=-\infty}^{\infty} x[n]z^{-n} < \infty $$

Now I'm confused I would say it doesn't exist because it has an infinite number of elements and its sum is infinite but then I check table for transform and transform from sigma function exist? It just jumps into infinity. Can someone tell me where I am making mistakes..


texts: why do directional derivatives work better with a gaussian filter?

texts:
I've been working on image segmentation and computing the directional derivatives of a grey-scaled image, with the objective of detecting contours and edges. I have realised that if I apply a gaussian filter before computing the directional derivatives, the edges get enhanced better. So my question is - why does that happen?

Here's the original image:

This is what the first order derivatives (dx and dy) look like (without gaussian filter):

here's the derivatives after applying gaussian filtering to the image, with sigma = 6:

*******Python code to compute the derivatives

with rio.open("vein_template.tif", 'r') as ds:
    RGB_arr = ds.read(masked=True)  # read all raster values
rgb_1 = np.rot90(np.transpose(RGB_arr, [1, 2, 0]))
#the rgb image contains 4 channels, so I'm changing their order, and making the image upright

%matplotlib
plt.imshow(rgb_1)

# Transform to grey-level img
gsi = 0.2989 * rgb_1[:, :, 0] + 0.5870 * rgb_1[:, :, 1] + 0.1140 * rgb_1[:, :, 2]
plt.imshow(gsi, cmap = "gist_gray")

#compute derivatives
dx = np.diff(gsi, axis = 0)
dy = np.diff(gsi, axis = 1)

f = plt.figure()
f.add_subplot(1, 2, 1)
plt.imshow(dx, cmap = "Greys")
plt.title("dx")
f.add_subplot(1, 2, 2)
plt.imshow(dy, cmap = "Greys")
plt.title("dy")


texts: How to determine the filter type of a discrete frequency spectrum

texts:
I have the discrete impulse response of a filter and i want to determine which kind of filter it is by using the DFT.

The impulse response is:

h0 = [0,1,2,1,0]

which I zero-pad for a better result

N0 = 10;
h1 = [zeros(1,N0),h0, zeros(1, N0)];

Then i created the N-Point DFT of it using FFT

H1 = fft(h1);

The picture shows the plot of abs(H1)

Now i know that k1 = 0 represents the steady component and at first I thought that this is the spectrum of a band-stop filter but now that I read a little bit more about the DFT I'm not quite sure because the DFT is always mirrored for real signals.

So my question is: Do I have to use only one half of the DFT spectrum to determine the filter characteristics? (Which would mean that this is a low-pass filter?)


texts: Deconvolution Using Complex Division in The Frequency Domain

texts:
Consider these two signals:

a = [1 1 0 0 0 0 0 0]
b = [1 0 1 0 0 0 0 0]

their convolution is

c = a * b = [1 1 1 1 0 0 0 0]

I am trying to obtain b by using complex division to divide the discrete Fourier transform of c by the discrete Fourier transform of a. I am aware that in general, there may not be a solution when attempting deconvolution in this way, due to division by zero issues etc. -- however, I want to understand why this doesn't work in this specific case.

The discrete Fourier transforms of all three signals are:

F(a) = [
    2                 (0)
    1.707 - 0.707 i
    1 - i
    0.293 - 0.707 i
    0                 (N/2)
    0.293 + 0.707 i
    1 + i
    1.707 + 0.707 i
]

F(b) = [
    2                 (0)
    1 - i
    0
    1 + i
    2                 (N/2)
    1 - i
    0
    1 + i
]

F(c) = [
    4                 (0)
    1 - 2.414 i
    0
    1 - 0.414 i
    0                 (N/2)
    1 + 0.414 i
    0
    1 + 2.414 i
]

The zero-th entry (marked (0)) is the sum of the number of 1 values in the signal, equal to $\Sigma_{i = 0}^7k_i$ for each kernel $k \, \epsilon \{a, b, c\}$; the value of the $N/2$-th entry (marked (N/2)) is equal to $\Sigma_{i = 0}^7k_i(-1)^i$.

The problem arises when F(c) is divided by F(a) at the $N/2$-th entry: this is 0/0, which is undefined, or is sometimes defined as evaluating to 1. However, the correct result for the complete division $F^{-1}(F(c) / F(b)) = F(a)$ can only be obtained if this instance of 0/0 in the $N/2$-th position evaluates to the value 2, not 1, since the $N/2$-th entry of F(b) is 2. Other than this value at the $N/2$-th position, all entries of F(b) can be recovered correctly by simply dividing entries in F(c) by the corresponding entry in F(a), using complex division.

Why is the N/2-th entry of the Fourier transform problematic in this way? Is there a robust way to derive what the correct result should be for the division of Fourier transforms at the N/2-th position? Or is this simply due to the fact that the value 2 is lost due to multiplication by 0 in the pointwise multiplication $F(a) . F(b)$?

Does this problem affect other components of the Fourier transform? Or is it only possible for this to happen at the N/2th position, which will be real-valued and may be zero?

Are there other gotchas to dividing signals in the Fourier domain in this way, or does this only happen in cases where there are zeroes in the Fourier domain of one of the factors?


texts: How does nearest neighbour, bilinear and cubic interpolation work in images?

texts:
More math is appreciated for each of the methods and references are appreciated. I have tried understanding from Wiki and matlab link but don't understand how the translation matrix is being used there.

What I seek is to be able to teach this to someone else. Thanks!


texts:
LTI system with impulse response:

$h(t) = e^t \sin(-3t) \text{u}(-t)$

I know the rule and formula , but I am lost in how to do this because of the absolute. Please help.

Sorry if I posted in the wrong section.

Thank you.


texts: Natural frequencies of a transfer function

texts:
I have been told in the university that the natural frequencies (also called $\textit{eigenfrequencies}$) are the poles of the transfer function, however, Matlab compute them as the modulus of the pole (information taken from here).

So the real question is, are the natural frequencies complex numbers or real numbers measured in rad/s?


texts: Interpolation algorithms in computer graphics

texts:
I'm looking for an image-processing algorithm that uses interpolation (polynomial, spline, doesn't matter). I have to implement such algorithm to process a bitmap for my labs. It can't be image resizing, it must be something else. Could you give me some proposals?
I've tried on wikipedia and google, but it's hard to filter algorithms that make use of interpolation

PS. Sorry for the tags, I can't create new ones and couldn't find existing ones that would be suitable for my question.


texts: Low pass filtering using MATLAB; passing frequency in Hz

texts:
I want to implement a HPF which will pass all frequencies above 150hz.

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

plot(Fs*t(1:50),x(1:50))
title('Signal')
xlabel('time (milliseconds)')

pause;

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
X = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
plot(f,2*abs(X(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of x(t)')
xlabel('Frequency (Hz)')
ylabel('|X(f)|')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Now let us see HP Filtering of this signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

d = fdesign.highpass('Fst,Fp,Ast,Ap,Fs', 150, 170  ,60,1, Fs);

Hd = design(d,'equiripple');
%fvtool(Hd);
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
HP_Filterd_Output = filter(Hd,x);

NFFT = 2^nextpow2(L); % Next power of 2 from length of x
HP_Filtered_Freq = fft(HP_Filterd_Output,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum of HPF.
plot(f,2*abs(HP_Filtered_Freq(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of HP_Filtered_Output')
xlabel('Frequency (Hz)')
ylabel('|HP_Filtered_Freq_Amplitude|')

I am getting following errors:

The 'Fst,Fp,Ast,Ap,Fs' enumerated value is invalid.

Error in ==> C:\Program
Files\MATLAB\R2010a\toolbox\signal\signal\@fdesign\@abstracttypewspecs\schema.p>set_specificationtype
at 65

Error in ==> fdesign.abstracttypewspecs.abstract_setspecs at 25
    set(this, 'SpecificationType', varargin{1});

Error in ==> fdesign.abstracttypewspecs.setspecs at 8
abstract_setspecs(this, varargin{:});

Error in ==> fdesign.highpass.highpass at 106
this.setspecs(varargin{:});

Error in ==> fftExample at 75
d = fdesign.highpass('Fst,Fp,Ast,Ap,Fs', 150, 170  ,60,1, Fs);

I guess I am making some mistake in

d = fdesign.highpass('Fst,Fp,Ast,Ap,Fs', 150, 170  ,60,1, Fs);`

where I want to enter every thing in Hz.


texts: What is difference between Phase shift keying (PSK) and Quadrature Amplitude Modulation (QAM) mathematical equation

texts:
I went through a literature and found the following expression for PSK signal. On this basis, what will be the equation for **QAM signal **. Also from where i would be able to read about this equation of PSK and QAM. Any help in this regard would be highly appreciated.
$$s\left[n \right] = \sqrt{P_s}\exp \left(j \frac{2\pi k}{M} \right)$$


texts: OCR (optical character recognition) for Indian languages - segmentation

texts:
Can any one help me in how to extract a table from the whole document image and segment it? I have done segmentation for the document image which has only text by using the projection profile. So please help me how to extract the table from the whole document image. I want only table to be extracted from the document. Please give me an idea how to do it. I want to continue my project using the projection profiles only.

i have done the segmentation manually for the document image which was consisting only text using projection profile manually i din't use the tables in it. now i want to extract the contents of the table in an image without using the formulas or any transforms can you help me.


texts: Least Mean Squares Algorithm confusion with Adaptive Line Enhancement

texts:
I'm just a bit confused about the least mean squares algorithm to separate wideband and narrowband in an adaptive filter for voice conversation. I'm interested in the narrowband part and I'm confused about the LMS equation as follows:

Is h being iterated twice here with j and n? Wouldn't that mean I have two different sets of h values?

I'm having trouble putting this into MATLAB.


texts: Does the Quality of the Pictures Taken by a USB 3.0 Camera Depend on Cable Length?

texts:
I've got a JAI GO-5000-USB camera and have tried several cable lengths, ranging from 40cm to 5m.

So far, I haven't observed many differences in the quality of the pictures I've taken. Since I want to avoid a later realization that the quality isn't that good after all, I'm hoping that someone with experience with USB 3.0 cameras and longer cable lengths could warn me in advance.

Now my question is:
Has anyone experienced a decreasing quality of taken pictures depending on longer cable lengths? When did the problems start and in which way did they appear? Is it possible to measure or visualize those effects?

I've tried to visualize them with histograms, but that didn't really show anything useful at all.


texts: Can time-invariance be determined from a given a transfer function?

texts:
I've the following function.
$$ G(z) = 2 + \frac{-1+5z^{-1}}{(1-0.5z^{-1})(1-z^{-1})}$$
Calculating it's inverse using $\mathcal Z$-Transform, I get the following function:
$$g[n] = 2\delta[n] + 8u[n] - 9(0.5)^{n}u[n]$$
where $u[n]$ is the unit step and $\delta[n]$ is the impulse.
I tried to do the usual approach to find if a system is time-invariant or not.
That is shifting
the input and compute the output $y_1$ and shifting the output $y_2$, compare both; if equal then
it's time-invariant, else it isn't time-invariant. However, I don't know what could be the input in this case.
So, I'm supposing that $g[n]$ is the impulse response of the system. But I couldn't find a way to tell if the system is time-invariant or not based on it's impulse response. Maybe I'm missing some
key idea here.


texts:
Can anybody explain the meaning of Doppler tolerance for the waveforms(Linear frequency modulated, Barker code, Polyphase codes) used for Pulse compression.

The Doppler frequency shift $f_{d}$ due to a moving target of velocity $v_{r}$ is given by

$f_{d} = 2*v_{r}/\lambda$ where $\lambda$ is the wavelength of the transmitted signal.

Polyphase codes are known to be having high Doppler tolerance compared to LFM waveforms.
Can anybody explain this fact in simple terms.


texts: What is the main difference between the acoustic and electromagnetic waves?

texts: For simulation purpose, what will be the difference between using acoustic and electromagnetic waves? Just the carrier frequency only??


texts:
In order to know if my signal is increasing or decreasing, I'm using the discrete derivative $y[n] = x[n] - x[n-1]$
or a smoothed version of it (for example Exponential Weight Moving Average of $y[n]$ that gives a "smooth" version of the derivative).
Now to make further analysis, I also want to consider the discrete acceleration:
$$z[n] = x[n] - 2 x[n-1] + x[n-2]$$
But as the signal is sometimes a bit noisy, I wanted to study a smoothed version of this acceleration. I thought about :

Standard Moving average of $z$ : $SMA(z[n], 10)$. Easy computation shows that $SMA(z[n], 10) = \frac{1}{10} (x[n]-x[n-1]-x[n-10]+x[n-11])$
It's not so good because it uses too few values of x[k]...

Exponential Moving average of $z$ : $EMA(z[n], 10)$. Computations shows that: $EMA(z[n], 10) \approx 0.09 (x[n] - x[n-1])$ and this is proportionnal to ... derivative !
Indeed, all the next coefficients are negligible $\ll 10^{-3}$.
This is strange: I expected the the Exponential Moving Average of the acceleration to be smoother but it turns out that it nearly uses only 2 values : $x[n]$ and $x[n-1]$. (See below the impulse reponse of EMA of acceleration).

Question:
How to get a "smoothed" version of discrete acceleration, as $A[n] = a_0 x[n] + a_1 x[n-1] + ... + a_k x[n-k]$ ?

Remark: Impulse reponse of EMA of discrete acceleration:


texts:
Rakeness optimization problem demonstrate that increases the rakeness between a , b  while leaving b random enough.

where e is the energy of the projection waveforms and  r is a randomness-enforcing design parameter.

How to achieve enough randomness in MATLAB?


texts: Reason for >8 channels in an audio codec

texts: I'm new to audio (and signal processing in general), and I've been learning about the AAC codec. I've read that it supports 1-48 synchronous channels in a single stream. My question is: what need is there for 48 channels in a single audio stream? I understand that if you have a 7.1 speaker setup you can send a separate channel to each speaker, but what is the need for 48 separate channels?


texts:
I am working on some experimental data which, at some point, need to be time-integrated and then high-pass filtered (to remove low frequency disturbancies introduced by integration and unwanted DC component).

The aim of my work is not related to filtering, but still I would like to analyze more in detail the filters I am using to give some justification (for example to motivate why I chosed to use a 4th order filter instead of a higher/lower one).

This is the filter I am using:

delta_t = 1.53846e-04;
Fs = 1/delta_t;
cut_F = 8;
Wn = cut_F/(Fs/2);
ftype = 'high';
[b,a] = butter(4,Wn,ftype);
filtered_signal = filtfilt(b,a,signal);

I already had a look here: https://stackoverflow.com/questions/5591278/high-pass-filtering-in-matlab to learn something about filters (I never had a course on signal processing) and I used

fvtool(b,a)

to see the impulse response, step response ecc. of the filter I have used.

The problem is that I do not know how to "read" these plots.

What do I have to look for?

How can I understand if a filter is good or not? (I do not have any specification about filter performances, I just know that the lowest frequency I can admit is 5 Hz)

What features of different filters are useful to be compared to motivate the choice?

For example, what info can I see plotting the round-off noise power spectrum of 2 different filters?(high pass Butterworth of different orders: 4 and 22)

EDIT

I tried decimating the signal to reduce the sampling frequency as suggested in the comments and this is what I get

The green one is the new filter's round off noise spectrum having reduced the sampling frequency (now = 325 Hz), the blue one is the older with sampling frequency = 6.5 kHz.

Things look a bit better but still not good enough, am I right?


texts: Upsampling question- Aliasing/imaging problem

texts:
I'm trying to study vibrations in a non-stationary shaft by upsampling the original signal.

The original sample rate is at 25.6kHz, and it is being upsampled to 36kHz. In the FFT of the upsampled signal I am getting freq. components above the Nyquist frequency of 12.8 kHz. I'm using linear interpolation when upsampling.

Is this imaging/aliasing? What are common techniques for dealing with such behaviour? I'm working in python/matlab.

Would scipy.signal.resample_poly.html be sufficient?


texts:
What is difference between edge detection and contour detection?

Apparently both techniques are used to detect boundries of close shapes(wheter oval shape or r rectangular shape)


texts:
I would like to implement a Rayleigh fading channel for an ASK Modulation. My question is concerning the receiver. How to properly add every tap at the receiver to keep the same length from the emitter.

For instance :
Emitter : binary frame : 1*N
Channel : six taps : M*N (where M is the number of taps)
Receiver : 1*N

Thanks.

Best regards.


texts: How to Calculate the number of DFTs taken in a STFT

texts: Given variables of window size (W), overlap (O), and length of the Signal (S) (in my case, length of a song) and I am trying to determine an equation that gives me the number of DFTs that would be taken in an STFT as the window slides given those variables above.  I'm not sure if those variables is enough to give me the equation or not, and if it is, I'm not sure how to derive that equation.


texts: autocorrelation and ISI histogram of oscillatory neuron

texts: I'm trying to learn about different kinds of neurons and I got to a question that I'm stuck on: Given an oscillatory neuron that shoots poissonically for 50 ms and rests for 100ms:how it's autocorrelation looks like and how does it's isi ( ISI= distance between spikes)histogram looks like?


texts: Optimal segment length for coherence estimation

texts:
I have a question about the computation of coherence between two signals with Welch's method.
The two signals are relatively short (e.g. 256 samples) and I would like as correct an estimate of shared frequencies between the signals as possible.
I am aware that I need to average over more than one segment of data for meaningful results and that overlapping segments will improve accuracy. Computational effort in the computation is not an issue, so I am using $noverlap = nseg-1$. Acquiring more samples is not possible.

The coherence is computed as $C_{xy} = \frac{|G_{xy}|^2}{G_{xx} G_{yy}}$ where $G_{xy}$ is the Cross-Spectral Density and $G_{xx}$ and $G_{yy}$ are the Power Spectral Densities of both signals, respectively.

My question is:
What is the optimal length of segments?
nseg=ndata (no averaging) leads to $C_{xy} = 1$ whereas a too small value of nseg leads to insufficient frequency resolution.
Is there any value for nseg for which an optimal compromise between sufficient averaging and frequency resolution can be found?

Both theoretical and best-practice suggestions would be appreciated!


texts:
Anyone know what references might have information about the history of any particular filter design?

Since seeing just the equation or implementation hinders the process by which the filter (such as the biquad filter) was actually discovered.


texts: Signal-to-Noise ratio of multivariate stochastic process from Correlation Matrix

texts:
I'm not in signal processing, I'm from an another discipline. I've derived a simple result which I presume must be well known in SP and I'd like to know whether there's a paper or textbook that has it that I can cite.
Suppose I have an n-dimensional stochastic process
$\mathbf{x}_{t}\;=\;\left[\begin{array}{cccc}
x_{t}^{\left(1\right)} & x_{t}^{\left(2\right)} & \cdots & x_{t}^{\left(n\right)}\end{array}\right]$
where $t$ denotes an observation. Each of the components $j$ is the combination of a signal $x_{t}^{*}$ and some random idiosyncratic noise:
$x_{t}^{\left(j\right)}\;=\;x_{t}^{*}b_{j}+\varepsilon_{t}^{\left(j\right)}\sqrt{1-b_{j}^{2}}$
where $\varepsilon^{(j)}_{t}$ is iid distributed and uncorrelated across the $n$ components.
For convenience, all these  variables are already standardized (mean 0 variance 1).
$b_{j}$ and therefore the signal-to-noise of each $j$ component ( $b_{j}/(1-b^2_{j})$ ) are not observed.
The simple result is: you observe the correlation matrix $\mathbf{R}$ of the $n$ components of $x_{t}$. Then a simple quadratic system gives you the noise-to-signal ratios of each $j$ component of $\mathbf{x}_{t}$ .
Anybody can tell me whether you've seen this result in a paper or textbook?
Thanks!


texts: Is there a difference measure that is insensitive to shifts?

texts:
I'm trying to figure out a difference measure that is relatively insensitive to time shifts.
I have tried DTW (dynamic time warping.)
This is the result:

The signals on the left are clearly more similar than the signals on the right yet due to time shifts I get larger normed difference for the left than for the right.
Is there any measure/technique I can use to overcome this? I will investigate the frequency domain but I wanted this community's opinion.


texts: Analytical Calculation of the Impulse Response of an Arbitrary Filter

texts:
This is a rather general question.  Suppose I have an arbitrary filter, $F$, defined in the Laplace domain;
$$
F=k\frac{\prod_i (s+z_i)}{\prod_j(s+p_j)}.
$$
How do I calculate the impulse response of this filter analytically?  As a specific simple example, consider the first order low pass filter
$$
F_0=\frac{1}{s+p}.
$$
To be clear, what I want to get from this transformation is a function which describes the delta function response from $t=0$ to $t\rightarrow\infty$ where the delta function is applied at $t=0$.


texts: What is a "Unit Shock" in an Impulse Response Function?

texts: Is a "one unit shock" in an impulse function of variable "temperature" a 1% increase or 1 more "unit" (1 degree)?


texts: Understanding multirate filtering fundamentals

texts:
I am having trouble grasping some of the fundamental concepts of multirate filtering. I see from various sources that the basic building blocks of a multirate filter are the dyadic analysis and synthesis blocks.

Question 1:

The analysis block structure looks like the following, where the wide band signal is split into lowpass and highpass bands, each with a cutoff of FS/4 (Nyquist/2). Each band is then decimated by a factor of 2.

How can you accurately represent the signal in the high frequency band when it contains frequency information above the Nyquist limit of the new decimated sample rate?
Question 2:

The analysis block structure looks like the following, where the sub-band signal is interpolated, re-filtered and then summed.

What is the purpose of the second filtering?


texts: Transform a Continuous Time Differential Equation of a Low Pass Filter to Discrete Time

texts:
Assuming I have a linear differential equation of first order low pass filter:
$$RC*y(t)'+y(t) = x(t) $$
And Assuming i know now that the input x(t) is being sampled at a rate of   $T_s$

how can now switch the equation to describe in the discrete time form?
I know that the answer is:
$$x[n] = y[n] + RC/T_s* (y[n]-y[n-1])$$

but i can't figure how to get there.
I know that there is a strong relation between DTFT and fourier transform but i couldn't figure how to use it.
Thanks


texts: What is the default carrier frequency of matlab's pskmod?

texts: I'm using the matlab function pskmod for modulating a bit string. What is the carrier frequency used by default by the function?


texts: Doubt on behaviour of frequency responce of a system

texts:
I am doing loop shaping using $H_{\infty }$ approch, and I have a plant with an unstable pole:

G = 10/((s+10)*(s-1));
K1 = 4*(s+1)/s;

where G is the plant, and K1 is the controller. I am trying to do a confrontation between the approch with $H_{\infty }$ and the one without, in particular I am placing a weight only for the sensitivity function, and designing a controller using the mixed sensitivity design.

What I get if I plot the frequency response of the system in the two cases is the following:

in which the yellow line is the frequency response of the system in which I applied the $H_{\infty }$ and the red line is the one in which I used G and K1 as they are, so without mixed sensitivity and weights.

What I don't understand is why the yellow line has this strange behaviour, while the red line decrease smoothly.

I see that the bandwidth of the yellow line is lower the the other one, so maybe this fact is due to the fact that I have a RHP pole which imposes a lower limitation on the bandwidth of the system. But, if I plot the sensitivity functions of these two systems, I have:

so the yellow line behaves better with respect to disturbance rejection with respect to the red line, as expected.

I don't understand how to interpret the behaviour of the yellow line of the complementary sensitivity function.

Can somebody help me?


texts:
I'm trying to do Analysis of EEG Signal for the Detection  of epileptic seizure so in the pre-processing of the EEG signal I use a low pass FIR filter

d1=designfilt('lowpassfir','PassbandFrequency',0.15,'StopbandFrequency',0.2,'PassbandRipple',3,'StopbandAttenuation',60,'DesignMethod','equiripple');
a = filtfilt(d1,deriv1);

I want to know if the values above for PassbandFrequency and StopbandFrequency are right.


texts: Trouble Generating a Sine Wave in C/Matlab from Basics

texts:
this is a crosspost from stackoverflow as I figured there may be more who can help me here...

I've been writing a basic C program to generate a sine wave and write it out to STDOUT for piping into a filter which I am also going to write but I'm having issues.

I am able to generate what sounds like a good clean sine wave, but the frequency is off.  When I listen to the output I am generating and compare to a true signal generators output, the frequency of mine is slightly too high.

The C code should generate raw 16 bit signed stereo wav @ 44100Hz for a frequency passed as an argument to the program, or 500Hz as default.

A matlab version of the code is here, (modded very slightly for 1 indexing in Matlab, WHY MATHWORKS WHY!?) minus the passing to STDOUT etc as I know that works fine

CHANNELS = 2; SAMPLING_RATE = 44100; NUM_SAMPLES = 512;

frequency = 1000;

%% BEGIN WHILE(1) LOOP HERE

output_buff = zeros(1,CHANNELS*NUM_SAMPLES);

for i = 1:2:CHANNELS*NUM_SAMPLES
    output_buff(i) = 30000 * sin(frequency * pi * (i-1)/NUM_SAMPLES);
    output_buff(i+1) = output_buff(i);
end

%% OUTPUT TO STDOUT AND REPEAT

I should add that this code runs (in the C version) inside a while true loop, generating a full output_buff of values, then pushing the buffer to STDOUT.

I have written some further testing code to see what's actually being generated as follows:

plot(1:CHANNELS*NUM_SAMPLES, output_buff)

output_buff = output_buff .* hanning(length(output_buff))';
Y = fft(output_buff);
Mag=abs(Y(1:length(output_buff)/2)).^2;
[a,b]=max(Mag);
% Result
SAMPLING_RATE*b/length(output_buff)

When I run this script, I can see at the end the frequency of the generated signal is actually 1.0767e+03Hz... Close but no cigar...

I've tried tweaking some of the parameters but I have no idea what's wrong or how to make the generated frequency more accurate.

The C code is as follows:

#include <sys/types.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>

#include <unistd.h>

#define M_PI 3.1415926535897932384626433832795L
#define NUM_CHANNELS 2
#define NUM_SAMPLES 1024
#define SAMPLE_RATE 44100

double frequency = 500.0;

int main(int argc, char *argv[])
{
    if (argc >= 2){
        frequency = atof(argv[1]);
    }

    int16_t samples[NUM_CHANNELS * NUM_SAMPLES];
    unsigned cbBuffer=sizeof(samples);

    while(1){
        for (int i = 0; i < cbBuffer; i+=NUM_CHANNELS){
            samples[i] = (int16_t) 30000 * sin(frequency * M_PI * (double)i /(double)(NUM_SAMPLES));
            samples[i+1] = samples[i];
        }

        int done=write(STDOUT_FILENO, samples, cbBuffer);
        if(done<0){
            fprintf(stderr, "%s : Write to stdout failed, error=%s.", argv[0], strerror(errno));
            exit(1);
        }else if(done!=cbBuffer){
            fprintf(stderr, "%s : Could not read requested number of bytes from stream.\n", argv[0]);
        }
    }

    return 0;
}

UPDATE:

I changed the C code as follows and I now get a nice smooth sine wave but the frequency is way too low

int main(int argc, char *argv[])
{
    if (argc >= 2){
        frequency = atof(argv[1]);
    }

    int16_t samples[NUM_CHANNELS * NUM_SAMPLES];
    unsigned cbBuffer=sizeof(samples);
    int counter = 0;

    while(1){
        for (int i = 0; i < NUM_CHANNELS * NUM_SAMPLES; i+=NUM_CHANNELS){
            samples[i] = (int16_t) 3000 * sin(frequency * M_PI * (double)counter /(double)(SAMPLE_RATE));
            samples[i+1] = samples[i];
            counter++;
        }

        int done=write(STDOUT_FILENO, samples, cbBuffer);
        if(done<0){
            fprintf(stderr, "%s : Write to stdout failed, error=%s.", argv[0], strerror(errno));
            exit(1);
        }else if(done!=cbBuffer){
            fprintf(stderr, "%s : Could not read requested number of bytes from stream.\n", argv[0]);
        }
    }

    return 0;
}


texts: A State space model for discrete Sine wave Using kalma filter

texts:
I'm looking to apply an optimal LQR filter to a discrete signal of the form
$$
x[n]=A\sin(\omega_0n+\phi)+v[n]
$$
The amplitude $A$ and the phase $\phi$ are unknown variables I want to estimate using the filter, and $v[n]$ is an uncorrelated noise signal of variance $\sigma_v^2$.

I don't know how to build a state model to generate this sine wave, and proceed from there.


texts:
How to calculate Hilbert transform, if it exists, of the signals like $u(t)$, $sgn(t)$.

What properties should a function satisfy for existence of Hilbert transform.

Absolute integrability of a signal $x(t)$ is necessary condition for HT to exist?


texts: Low pass filter implementation with programmable cut-off frequency for inclinometer

texts: My application is an inclinometer. I would like to design a low pass filter with programmable cut-off frequency ranging from 0.1Hz to 20Hz using CMSIS DSP library provided by ARM. The filter may be critically damped and butterworth with an order of 1, 4 and 8. The sampling rate of the sensor can be configured between 100Hz to 1600Hz. How can I implement it? Thanks in advance.


texts: How to generate Gray codes of n bits that have at max k ones?

texts:
I understand Gray codes have one bit change between them and how to convert between binary and gray. But given a bit length of n, I want to generate a series of Gray codes (all possible) that have one bit change and have maximum k ones.
Example: Given n = 3 and k = 2
000 001 011 010 110 100 101


texts:
I am trying to plot a spectrogram using Matlab.
Basically what I need is to have:
-time in the X axis
-Velocity in Y axis
-Intensity is the current (in pico Ampere)

Basically, I want to create a simliar plot as the one attached here.

Time is an array of 9000 elements
Velocity is a matrix of 9000x128 elements
Current is a matrix of 9000x1280 elements

I have a thermal velocity which I can also plot which is an array of 9000 elements.

Can anyone give sugggestion on how to plot it?
I want to crate similar to the one below:


texts:
I have a question regarding the fourier transform. If I plot the frequency against the fourier transform for a periodic signal and I get a peak, What is the physics behind it? I want to know the physical meaning of it.

Why we get a peak and what does it signify?

Thanks


texts: Discrete Fourier Transform: What is the DC Term really?

texts:
I am currently toying around with the Discrete Fourier Transform (DFT) in Matlab to extract features from images. I like to fully understand the concepts that I use. I have read several explanations, such as this, but so far, none really explained the meaning of the "DC term". All I know is that the k'the term of the DFT can be written as:

where
 is the twiddle factor.

That means that the first term (the DC term), , is an amplitude without frequency.

Can someone explain why is it called the DC term? What is it's relation to "Direct Current"? And what is the relevance of the DC term? When is it useful, and for what?


texts: Significance of modular arithmetic in DFT?

texts:
In what ways does modular arithmetic plays a part in DFT?
Why is it a so integral part of DFT?


texts: The magnitude of taking discrete cosine transform of an image two times is similar to the original

texts:
Consider the following:

im = double(imread('lena.bmp'));
subplot(1,3,1), imshow(im,[]), title('original');
t1 = dct2(im);
subplot(1,3,2), imshow(log(abs(t1)+1),[]), title('DCT transform');
t2 = dct2(t1);
subplot(1,3,3), imshow(t2,[]), title('DCT(DCT) transform');

the output is shown below:

Could anyone explain why the result of the second order DCT is similar to the original image?

Thanks in advance!

P.S. The same thing happens with every 2N'th DCT.


texts: Linear and Circular Convolution in Fourier Domain (DFT)

texts:
Suppose we have two vectors A and B of length 100 and 80 obtained as a function of time. If we wish to perform convolution of the two vectors in the Fourier domain, we need to multiply the Fourier transforms of A and B. For that purpose, the length of B must be the same as A. The output would be of the same length as A or B after zero padding.

Where should we do zero padding in B, should this be at the beginning of B or at the end of B?
I was reading that convolution achieved via FFT is essentially a circular one. If our desire is to get a linear convolution, how do we ensure that the output is a linear convolution, i.e., how many end points should be rejected from both ends after inverse FFT?


texts: Is it correct to subtract two signals acquired over different times and trials to remove the common signal in them?

texts:
I have two signals that represent the response of a neuron under two different conditions.

Signal 1 (S1): response to Stimulus A
Signal 2 (S2): response to Stimulus A+B

The response to stimulus A is common in both signals but I am interested in the response to stimulus B only. Because the signals are noisy, they have been acquired multiple times. S1 is acquired 60 times and S2 is acquired only 4 times.

Will averaging S1 over these 60 acquisitions and S2 over 4 (Figure 1) and then subtracting the averaged signals (S2_avg - S1_avg) give me a better representation of response to Stimulus B or should both the signals be averaged across equal acquisitions (Figure 2)?

Averaging S1 60 times really smooths the response to Stimulus A but compared to this, response to Stimulus A is not smooth in S2 as it is only averaged 4 times. So in this way, what I understand is  by averaging S1 many times might actually be a disadvantage.

  Additional info following answers: my underlying assumption is that
  the neuron behaves like an LTI system. My concern is that averaging S1
  many times smoothes it as compared to response to stimulus A in S2. As
  such, subtracting them would not really remove the common signal in
  them.


texts:
I am trying to understand how to use nufft from the Matlab doc.
My goal is to compute a FFT of an image (2D) which has missing points (not sampled). I have a list x and y of coordinates of the points that ARE sampled. These lists are for the x and y coordinates, like, say for a 64x64 image: x=[1,2,3,10,11,12...],y=[1,2,8,32,32,33,...]
t = [1:10 11:2:29]';
x = t;
y = t';
z = reshape(t,[1 1 20]);
X = cos(2*pi*0.01*x) + sin(2*pi*0.02*y) + cos(2*pi*0.03*z);
Y = nufftn(X,{t,t,t});

From this example I tought I could call the nufft() function like so:
 x=imread('cameraman.tif');
 xx=[1:2:256];yy=[1:2:256];
 xnufft=nufftn(x,{xx,yy});

However Matlab outputs an error message:
**Error using nufftn (line 107)
The total number of sample points must match the number of elements in the first input argument.**

I do not understand this message since by definition if it is "nonuniform" I will have less sample points in each dimension of the input ?!? I do
I do see that in the example, their X has the same number of values in each dimension but the problem is that I want to give as input an image with "holes", so unless I can put in a special value for the "holes" I don't see how to input the image unless I vectorize it and remove the "missing" samples but the nufftn is supposed to be designed to nD data.
I give as input an image img with zeros at the points that are not sampled, but I do not see how I could give an image with actual "holes" since it would result in a matrix with different number of columns in each row which is impossible.


texts:
How is it possible to go from time domain amplitude values in floating point to frequencies?

I know the DFT is used, but what exactly is happening for it to be able to extract frequencies from amplitude values?

Is amplitude between -1 and 1 mapped from pressure levels? How are different frequencies detected in a block of amplitude values? What if two files have the exact same amplitude values but completely different frequencies?


texts: Kalman Filtering with Unknown State Transition Matrix

texts: I'm currently studying the use of Kalman filters for estimating linear systems. My current State Transition Matrix (STM) is the identity since so far I've been dealing with non time-varying systems. However my next step is to extend to time-varying systems and in this case there is no information on how the system varies and so the STM is unknown. And as of my testing so far a wrong STM causes a ton of errors on the Kalman Filter estimations. So my question is if there is a method for circumventing this, be it by updating the STM on the fly or some kind of way of updating the Process Noise Covariance Matrix in order to account for that?


texts: Discrete Convolution of two piecewise sequences having this specific form

texts:
Assume I have the following two sequences :
$$
x[n]=\begin{cases}
\alpha&\text{if $a\leq n\leq b$}\\ \\ \tag{1}
0&\text{if otherwise}
\end{cases}
\qquad
\text{and}
\qquad
h[n]=\begin{cases}
\beta&\text{if $c\leq n\leq d$}\\ \\
0&\text{if otherwise}
\end{cases}
$$
I was wondering in this case if the convolution :
$$
(x*h)[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k] \tag{2}
$$
has a specific form applicable for these two piecewise sequences.
I have seen on the internet something like :
$$
(x * h)[n]=\sum_{k=\max (a, n-d)}^{\min (b, n-c)} x[k] h[n-k] . \tag{3}
$$
However, I can not confirm if this is true. I hope someone can provide a proper proof and thank you very much.


texts: Implementation of gain scheduling for a PI controller

texts:
I'm designing a PI controller for a boost converter. I was able to find a set of gains that fit for all situations, but I get a slow response at full-load.

I was thinking instead of using gain scheduling, 1 set of gains for little to no load and 1 set of gains when load is significant (say more than 25% of the maximum load).

Now my question. What are the gotchas about gain scheduling? In books they don't say much about implementation details.

I was thinking that in order for gain scheduling to work, I need to apply the Ki gain of the integrator AFTER the integration and not before.  That way I canchange seamlessly from one set of gain to the other.

I.e :

Use this =>
integrator = integrator + error;
output = Kiintegrator + Kperror;

Instead of
integrator = integrator + Kierror;
output = integrator + Kperror;

Am I on the right path?

Regards

Benjamin


texts: How to draw a vector between light and dark spots on an image

texts:
I'm very new to Python and I need to draw vectors between dark and light spots on an image. Ive managed to convert it to greyscale but I'm struggling with identifying the dark and light spots. If I can identify the co-ordinates of the centre of the bright spots then I should be able to draw a vector on. The aim is to get it to work on an image with one pair of spots then extend this to multiple pairs.

Thanks to anyone who helps.


texts: How can I implement a triangular filter in MatLab, given it's Z-transform representation?

texts:
I have to implement an anti-aliasing filter for a certain processing step in MatLab. While searching in literature for inspiration, I came accross a paper in which the authors wrote the following:

We implement the lowpass filters as triangular smoothing. The Z-transform representation of an arbitrary N-point filter of this type is:
$$G(z)=\frac{-z^{-(k+1)}+2-z^{k+1}}{\alpha(1-z)(1-z^{-1})}$$
where the triangle length $N=2k+1$, and the filter scaling coefficient $\alpha=(k+1)^2$.

I was thinking about how to implement this filter, so following the approach taken in this DSP SE post, I rewrote the expression for the filter as follows:
$$G(z)=\frac{-z^{-(k+1)}+2-z^{k+1}}{\alpha(1-z)(1-z^{-1})}=\frac{Y(z)}{X(z)}$$
$$X(z)\Bigl(-z^{-(k+1)}+2-z^{k+1}\Bigr) = Y(z)\Bigl(\alpha(1-z)(1-z^{-1})\Bigr)$$
let's assume for now we use a filter of length $N=3$, which gives us $k=1$, we get:
\begin{align}
X(z)\Bigl(-z^{-2}+2-z^{2}\Bigr) & = Y(z)\Bigl(4(1-z)(1-z^{-1})\Bigr)\\
X(z)\Bigl(-z^{-2}+2-z^{2}\Bigr) & = Y(z)\Bigl(4(2-z-z^{-1})\Bigr)\\
X(z)\Bigl(-z^{-2}+2-z^{2}\Bigr) & = Y(z)\Bigl(8-4z-4z^{-1})\Bigr)\\
-x[n-2]+2x[n]-x[n+2] &= 8y[n]-4y[n+1]-4y[n-1]\\
\end{align}
which I finally rewrite as:
$$y[n+1]=\Bigl(8y[n]-4y[n-1]-2x[n]+x[n-2]+x[n+2]\Bigr)/4$$
I implemented this result in matlab as follows:
n = 1:1000;
x = sin(0.10*n) + (-0.5+rand(1,1000));
out=zeros(size(x));

for n_idx = 3:998
    out(n_idx+1) = (8*out(n_idx) - 4*out(n_idx-1) - 2*x(n_idx) + x(n_idx-2) + x(n_idx+2) )/4;
end

subplot(2,1,1); plot(n,x);
subplot(2,1,2); plot(n,out);

but unfortunately, the output looks absolutely nothing like a lowpass filtered version of the input. The authors do not go in detail about how they implement this filter, the only clue they give is this line:

we note that the denominator of $g(z)$ represents a causal integration of the data with $1/(1-z)$, followed by an anticausal integration with $1/(1-z^{-1})$.

which I can't make sense of.
My questions are therefore:
-1 Is my implementation of the filter $G(z)$ correct?
-2 Would anyone know what the authors mean in their comment on (anti)causal integration?


texts: Multichannel input/output evaluation board

texts:
I'm working on a project at school and the first step is to choose a dsp/evaluation board but I'm having trouble. The project is to create a system that can be used for a variety of audio applications (such as sound localization or digital effects processing). Because of this the evaluation board should have at least 3 inputs and 1 output. Some other requirements at 16 bits/sample and at least a 16KHz sampling rate and the ability to operate in real time. The final result should be a setup where to do a new audio project, one would only have to do simple coding (preferably C, or maybe python). Thanks for any recommendations or general help, feel free to ask more questions about the project.

Ryan


texts:
I am getting the wrong output with the following simple implementation of Python fft.  Can someone please explain what I need to fix?

import numpy as np
import matplotlib.pyplot as plt

def tone(fs, frequency, dur=1, amp=50):
    '''
    test tone
    '''
    nsamples = dur * fs
    t = np.linspace(0, dur, nsamples, endpoint=False)

    tone = amp*np.sin(frequency*2*np.pi*t)
    return tone

fs = 16000
dur = 0.5
freq = 1500

stim = tone(fs, freq, dur=dur)
time = np.linspace(0, dur, num=fs*dur)

t = np.linspace(0, fs, num=fs*dur)
fdata = np.fft.fft(stim)

hwp = int(np.shape(fdata)[0]/2)
tplt = t[:hwp]
fplt = fdata[:hwp]

plt.figure(num=1)
plt.plot(time, stim, 'b-')
plt.show()

plt.figure(num=2)
plt.plot(tplt, np.abs(fplt.real), 'b-')
plt.xscale('log')
plt.show()

Removing

, endpoint = False

will fix the implementation though I do not understand why.
Perhaps someone can explain(?)


texts:
I am experimenting with the OpenCV-library in python. Let's say I read in an image such as:

import cv2
import matplotlib.pyplot as plt
import numpy as np

img = cv2.imread('imgs/watermark_no_copy.jpg')

This watermark_no_copy image should be part of some python-library.
What I want to do now is to extract just the part of the image that is not white. I have come across the following code to do the job:

 img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
 mask_inv = cv2.bitwise_not(img_gray)
 img_extracted = cv2.bitwise_and(img, img, mask=mask_inv)

I don't fully understand what's going on tho. I understand that we convert the image into a Grayscale at first. A grayscale-image consists of just one color channel (only black and white with different intensity). In my case the background of the image is white and should, therefore, be assigned intensity 255. The colored part of the image, however, is red. It will not be assigned intensity 0.

Now, how does bitwise_not work? I understand that this should flip any bit 1 to 0 and any bit 0 to 1. However, in that grayscale-image, we don't have just 0 and 1. We have different intensity values. How exactly, does this work? Same question can be applied to the bitwise_and in the last line. bitwise_and of an image with itself, should return just that image, shouldn't it?


texts: Is the fundamental basis function of the DFT not an intregral number of cycles?

texts:
When the DFT is defined and described I typically see the basis functions described as an integral number of sine and cosine waves.

But the basis functions are $\sin(2\pi k (n/N))$ or $\cos(2\pi k(n/N))$ where $k$ is the harmonic, $N$ is the number of samples, $n$ is the time domain sample number.

So the right hand boundary of this argument is  $2\pi k(N-1)/N$, which seems to be always a bit short of a full cycle.

Could someone please confirm or explain my error?

The picture below shows a sine basis function, apparently ending short of a full cycle.


texts: DSP or signal/image/data processing jokes

texts:
Some other StackExchange/StackOverflow sites are entertained with a certain level of humor or fun.  What is your favorite “data analysis” cartoon? is particularly worth  mentioning (IMHO)  this xkcd cartoon on causality and correlation (DSP folks know what I mean):

So Is humor allowed in answers? I hope so, se for instance:

StackOverflow: What is your best programmer joke?
MathOverflow: Do good maths jokes exist
SE.Stats (Cross Validated): Statistics Jokes
SE.Stats (Cross Validated): What is your favorite “data analysis” cartoon?
SE.Maths Questions about math jokes

Humor in signal processing or image processing does exist, albeit shallow, and narrow-spread. I have encountered it first (I had zero humor before) at the ICASSP 2000 DSP humor exhibition in Istanbul, Turkey. The image below comes from Humor in DSP by EURASIP:

So I believe that we deserve some humor and joke tags and actual lightness.
 Go ahead!


texts: How to choose the appropriate number of samples for discrete wavelet transform denoising?

texts:
I want to make sure I am understanding the filter bank frequencies correctly in DWT decomposition.
If my sampling frequency is 30kHz, then would I need to input data with 30k samples which later will be decomposed to:
(level 0 - 0-15k Hz)
(level 1 - 0-7.5k & 7.5k-15kHz)
(level 2 - 0-3.75k, 3.75k-7.5k & 7.5k-11.25k, 11.25k-15kHz)
and so on...
So if I need to zero power-line interferences (60Hz) but my signals of interest are in the low frequency range then I would need a min of 7 levels to get to 0-117 Hz and that still does not provide very high resolution in frequency (can eliminate other signals of interest within the range)
So my question is, does the number of samples for input as data to be decomposed have to be related to the sampling frequency or it can be any number below the sampling frequency to obtain higher frequency resolution?


texts: Offset after Butterworth bandpass filtering

texts:
guess I'm stuck with this filtering problem and I'm thinking it has something to do with my raw data but I'm not so sure. So the problem is the following, please refer to the figures below, basically, all I'm trying to do is implement a bandpass Butterworth filter to filter the raw data trace (green). The raw data consists of 2501 points and was sampled at 20 kHz (100-ms). I'd like to bandpass filter this trace between 3 Hz and 170 Hz using the Matlab function 'butter' followed by either 'filter' (top figure) or 'filtfilt' (bottom figure), preferably the latter.
My Matlab code is as follows:
A = Average(1,:); % raw data
fs_Hz = 20000; % sampling rate (Hz)
order = 3;
fcutlow = 3; % Hz
fcuthigh = 170; % Hz
[b,a]  = butter(order,[fcutlow,fcuthigh]/(fs_Hz/2), 'bandpass');
x = filtfilt(b,a,A);
plot(x);
hold on
plot(A,'color','g');
legend('filtered data','raw data');

As you might guess the problem is the offset that I observe after both filtering functions. These 2501 point traces actually stem from a very long data-trace, should I possibly filter this long trace prior to extracting these smaller traces? Is the reason for the offset the fluctuating pattern of the trace? I can't detrend the data as I need the y-axis information for the subsequent analysis on this data. Should I be using a completely different type of filtering method; FIR instead of IIR??
Any advice on my problem would be greatly appreciated.
IljaPilja


texts: Simulating two signals with a desired coherence

texts:
I'm trying to simulate two frequency-domain signals with a desired coherence, and I'd like to check if my idea is right.

To simulate two time-domain signals with a desired correlation, we can use the following method:

Suppose $\Sigma$ is the desired correlation. Let $X \sim N(0,1)$.

Then if we take

$X_1 = \Sigma^{1/2} X$,

$X_1$ will have the desired correlation:

$\mathbb{E} [X_1 X_1^{\dagger}] = \mathbb{E}[\Sigma^{1/2} X (\Sigma^{1/2} X)^{\dagger}] \\
= \mathbb{E}[\Sigma^{1/2} X X^{\dagger} \Sigma^{1/2 \dagger}] \\
= \Sigma^{1/2} \mathbb{E} [X X^{\dagger}] \Sigma^{1/2} \\
= \Sigma $

I'm trying to extend this to the frequency domain, after applying a Fourier transform to $X$. Let $X^F$ be the Fourier transform of $X$, ie. $X^F = WX$ where $W$ is the DFT matrix. Let the desired coherence be $S$. To use an idea similar to that used to simulate a desired correlation, I choose

$X_1^F = S^{1/2} S_{XX}^{-1/2} X^F$

where $S_{XX}$ is the cross-spectrum of $X^F$. Then

$\mathbb{E}[X_1^F X_1^{F\dagger}] = S^{1/2} S_{XX}^{-1/2} \mathbb{E}[X^F X^{F\dagger}] S_{XX}^{-1/2 \dagger} S^{1/2 \dagger} \\
= S^{1/2} S_{XX}^{-1/2} S_{XX} S_{XX}^{-1/2} S^{1/2} \\
= S$

I wrote a Python program to try this out but I didn't get a coherence close to my desired $S$, so I'm hoping to check whether there are any conceptual errors (otherwise it might be a numerical issue?). In particular, I'm not sure if $S_{XX}^{-1/2}=S_{XX}^{-1/2 \dagger}$ in general.


texts: Applying fourier transform twice (DSP course)

texts:
I stumbled upon a question in a DSP course (coursera) which I don't understand, shown below is a screenshot of the question and answer.

The part which I don't understand is circled. Why is it equal to 0? For multiples of N I understand that each element in the summation is 1 because the exponent is 0. But why does the summation equal 0 when (i+n) is not a multiple of N?


texts:
I'm very new to DSP hardware, so I am taking it upon myself to learn it with regard to Audio and Acoustic Signal Processing. I've learned quite a bit of theory in a graduate DSP course, but absolutely no practical implementation was introduced.

I've stumbled upon a site of miniDSP kits as well as a Texas Instruments TMS320C6748 Dev kit. What hardware should I purchase for learning audio DSP? And, if any, are there any good books that highlight the relative importance of different DSP hardware attributes? I would love to learn as much as I can about developing DSP hardware.


texts:
I have a simple question

How could one implement the $z$ or $z^2$ in simulink? The $z^{-2}$ could be implemented by Variable Integer Delay block, but what about $z^2$?


texts: Does flat surface detection and orientation require special hardware to compute?

texts:
Modern augmented reality platforms such as Google's ARCore and Apple's ARKit seem to only operate on mobile devices, I'm guessing, because their underlying algorithms require specialized hardware that is typically available on these devices (accelerometers, gyroscopes, etc.). Since this hardware isn't available on laptops/desktops, I'm guessing these libraries would never be able to work off of mobile platforms.

Having said that I'm wondering if the following AR/image processing "capabilities" can be achieved via algorithms that do not require such specialized hardware:

Flat surface detection inside an image or video; and
Orientation/angle of the camera view inside an image or video

Meaning, if I am given an image or video (set of sequential images), can I detect flat surfaces in that media without specialized hardware (as in, on a Linux box running on server or PC)? Same goes for orientation/angle determination (meaning, determining the (x,y,z) coordinates of the camera within the "domain" of the image/video and determining which angles its pointed towards)?


texts: Prewitt operator and central difference?

texts:
When taking the gradient of an image one can do $I * \begin{bmatrix}+1 \ 0 \ -1\end{bmatrix}$ (x direction).

As I understand this is basically applying central differences in the x direction, $\delta_h[f](x) = f(x+\tfrac12h)-f(x-\tfrac12h)$. But since the gradient is a vector of the derivatives and $f'(x) \approx \frac{\delta_h[f](x)}{h}$ shouldn't it be $I * \frac{1}{2} \begin{bmatrix}+1 \ 0 \ -1\end{bmatrix}$?

And what about the the Prewitt operator $\begin{bmatrix}+1 \ 0 \ -1 \\+ 1 \ 0 \ -1 \\ +1 \ 0 \ -1\end{bmatrix}$ (x direction)? As I understand there's some relation to finite differences as well, but how does having a 3x3 kernel affect this relation?


texts: Measuring sharpness/contrast for autofocus

texts:
For the implementation of an auto-focus for the raspberry cam I need some measure of image sharpness. I used the following sample picture:

The first approach was to compute the shannon entropy of the whole image over its histogram. Then I artificially blurred the image to see the change in entropy. Expectation was that the entropy monotonically increasing with increasing blur. But in practice this was not the case. The increase in entropy was very little, and at some point it starts decreasing again for some reason.
The fact, that the entropy increase is weak lies probably in the nature of the image: The majority of the image is already uniform. Blurring out the vessel does not change too much on the entropy.

If you are interested in the code, here it is:

Created on Tue Jun 28 23:54:58 2016

@author: asus
"""
import cv2
import cv2.cv
import numpy as np
import time

cv2.namedWindow('output')
src = cv2.imread("3.jpeg",1)
img = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)
#img = cv2.medianBlur(img, 5)
for i in range(1,300):
    entropy=0
    temp = cv2.blur(img, (i,i))
    hist = np.bincount(temp.ravel(),minlength=256)+1
    hist=hist/float(np.sum(hist))
    entropy=np.sum(-hist*np.log2(hist))
    print entropy

#    circles = cv2.HoughCircles(temp, cv2.cv.CV_HOUGH_GRADIENT,1,10)
#
#    if circles is not None:
#
#        circles = np.round(circles[0, :]).astype("int")
#        for (x, y, r) in circles:
#            cv2.circle(temp, (x, y), r, (0, 255, 0), 4)
#            while not cv2.waitKey(1)&0xFF==27:
#                    cv2.imshow('output',temp)

while not cv2.waitKey(1)&0xFF==27:
    time.sleep(0.1)
    cv2.imshow('output',temp)

cv2.destroyAllWindows()

As you may have noticed, I also tried to apply the histogram method only to few sections of the image. In this case I searched for circles via Houghtransform, and then only histogrammed the rectangles around the circles. But the circle detction itself is not really reliable for blurred images.

How can I enhance this method?
What other options do I have to obtain a measure of sharpness?

EDIT: After overthinking this, the decrease of entropy after a certain degree of blur actually makes sense: Tf the image is too blurred to strongly, then the image will be more monochrome, which means a peak in the histogram--> entropy converging to zero.

EDIT:

In response to A_A's comment I tried the following idea:
devide the image into 20 segments. blur the image at several different strengths (like going through all camera foci). For each of the 20 segments observe, how the entropy changes as a function of blur strength (or as a function of camera focus). Those segments, which do not contain too much of the objects will not change too much with increasing blur. Their entropies even decrease. Those with object information will have increasing entropy. Select those segments, and apply the entropy method on them all at once.

Results are little better. I will try that also with standard deviation instead of entropy, and then select the relevant segments with a threshold.


texts: Improving accuracy of Geortzel using simultaneous frequencies

texts:
I have written software to detect the transmission of frequency sequences.
I'm using Geortzel's Algorithm to filter desired frequencies from a data array. The data array is generated from audio capture. I am transmitting a sequence of bits encoded in two frequencies, F-mark and F-space.

Currently, my signal-to-noise ratio is not satisfactory and even with scaling and thresholding of the signals, I am getting too many false positives.

I've had the idea to use two independent frequencies to represent each the mark and the space (4 frequencies in total) and only register a valid signal if both are detected simultaneously. This seems like it would only approximately double the computation time which is acceptable.

Would this improve improve the accuracy of my detection, or would it have no effect on my signal-to-noise ratio? Would it also increase the likelihood of false negatives and missed detections?


texts: improve delay detection between two signals

texts: Hi there I'm trying to measure time delay of one percussive event recorded by two contact mikes on a glass surface. I decided to check the cross correlation (via FFT) and detect the index of the max in the function. The method is really jittery as the medium affects the transfer function of the signal. I wonder if there is a better method than cross correlation or some signal cleaning that can help before the calculation, or in the cross-correlation function. Also there might be a way to reconstruct the distortion transfer through the medium (as in phase cancellation)..any help?


texts:
I'm just beginning to learn about Fourier series and I'm trying to figure out how to find the Fourier series coefficients for $$x(t) = e^{j100\pi t}$$

I know that $$x(t) = \sum_{-\infty}^{\infty} a_{k} e^{jk(2\pi/T)t}$$

How do I go about finding the coefficient(s)? I know that I can get T = 1/50, but beyond that I don't even know where to begin. I think I'm supposed to be able to do this just by looking at it without having to solve the integral equation for $a_k$, but I don't know what I'm supposed to be looking for.


texts: Frequency response Test (Bode Diagramm) on system controlled with a PLC

texts:
i have to perform frequency response test on a pump system with servo Motor controlled by a PLC von ALlen Bradley (Logix 5156). I not sur how to proceed and will like to know if someone did such work before ? what is the best approach ? did i need a FRequency REsponse ANalyzer or it is possible to use the PLC to send sine wave to the system ? Or does a signal generator with an Oscilloscop will do the work ?
Thanks guys for your reply!


texts: Obtaining Correct (Least Squares Sense) Affine Transform Parameters Between Two Images

texts:
I have two images that I want to compute the affine motion model parameters. The model that I use is $$x' = a_1x+a_2y+a_3$$ $$y' = a_4x+a_5y+a_6$$
To calculate those 6 parameters, I picked 6 points (overdetermined system) between two images, and calculated those parameters by using MATLAB (A = X\X_primes). However, when I apply the transformation, even the points that I picked (X) do not translate into the specified locations (X_prime). I inferred that my least square solution should be wrong, and I am looking for how to make it better, i.e. what kind of points should I choose on the image to obtain better parameters for affine motion model?


texts: Implementing blind deconvolution in MATLAB

texts:
I want to implement blind deconvolution for the signal $ r(n) = h(n) \star s(n) + a(n)$ in MATLAB where

$r(n)$ is the recorded speech
$h(n)$ is impulse response of room acoustics
$s(n)$ is desired speech signal
$a(n)$ is noise from microphone

I understand in order to find the desired speech signal, $s(n)$, I need to perform some sort of blind deconvolution. I've been sent all over the place and was looking for some thoughts on how this could be approached.

Note: This is a university bonus project and I am in no way asking for an entire solution. I'm just looking for thoughts and maybe some hints about sites I can consult.


texts: The expectation in power spectral density

texts:
I'm a bit confused with the definition of the power spectral density (PSD). From Wiki https://en.wikipedia.org/wiki/Spectral_density , I found the definition is:

$$
S_{xx}(\omega) = \lim_{T\rightarrow \infty}\mathbb{E}[|x(\omega)|^2],
$$
where $x(\omega)$ is the Fourier transform of the process $x(t)$.

I'm really confuses with that expectation $\mathbb{E}$ in the above equation. The expectation is taken w.r.t. $x(t)$, but $|x(\omega)|^2$ is a function of angular $\omega$. There is no $x(t)$ in the integrand at all, thus
$$
\mathbb{E}[|x(\omega)|^2] = |x(\omega)|^2.
$$
What is the point of this expectation?

However, in the other hand,
$$
\mathbb{E} \left[ \left  | x(\omega) \right |^2 \right] = \mathbb{E} \left[ \frac{1}{T} \int_0^T x^*(t) e^{i\omega t}\, dt \int_0^T x(t') e^{-i\omega t'}\, dt' \right] = \frac{1}{T} \int_0^T \int_0^T \mathbb{E}\left[x^*(t) x(t')\right] e^{i\omega (t-t')}\, dt\, dt' \neq |x(\omega)|^2
$$
the PSD is the Fourier transform of the cross-covariance (auto-correlation) of the process. This exectation is indeed needed.

Where am I wrong? Can't those integrals switch orders?


texts:
I have dt signal $x[n]={[6.29, 8.11,-7.46,8.26,2.64,-8.04,-4.43,0.93,-9.29]}$

And I need to give the function value of:

1) sum of $x[n] =  \sum\limits_{k=0}^{N-1} X[k] \, e^{-j 3 \pi k/5}$ from k=0 to 9

How can I find that


texts: Does SWT/ISWT require intermediate approximation coefficients to represent/reconstruct the original signal?

texts:
Taking into account that Stationary Wavelet Transform (Algorithme à trous) is not an orthogonal transform do we need intermediate approximation coefficients for signal decomposition?

For example, in maximally decimated DWT we have the following decomposition tree:

And according to the conclusions from the MRA (multi resolution analysis) to fully represent (and reconstruct) the signal it is suffice to have only the last approximation subband plus all of the difference subbands.

Does this also hold true for SWT?


texts: How to get more accurate results with OpenCV

texts:
I would like to ask about this function  in OpenCV:

cvTermCriteria( CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, .3 )

How do I get more accurate results - do I need to increase the value of epsilon or decrease it?


texts:
I have some FFT data, 257 dimensions, every 10 ms, with 121 frames, i.e. 1.21 secs. I think the first dimension is probably something else and the remaining are the FFT coefficients, I guess.

From there, I want to calculate back some PCM signal for 44.1 Hz so I can play the sound. I asked the same question in a more mathematical way here but maybe here is a better place because I actually want to implement this.

How would I do that? Maybe I need some more information (which I have to find out somehow) - which? Maybe these missing information can be intelligently guessed somehow? What is a likely interpretation of those 257 dimensions?

This question is about the theory. The implementation is trivial I guess. But a concrete example in some language would be nice to help understanding the theory. Maybe C++ with FFTW? I skipped through the FFTW docs but I fail to understand all the terminology and some background, e.g. here. Why is it from complex to real or the other way, I only want real to real. What are those REDFT? What's a DCT, DFT, DST? FFTW_HC2R?

(I asked mostly the same question on SO here, but that was intended to be more about the practical side. Here on DSP I want to ask about the theory to get some understanding.)


texts: Estimating confidence intervals for the mean value

texts:
Say I conducted a numerical experiment -a time-dependent fluid simulation- and extracted a time-series of some variable. It looks like this:

As a result of the simulation, I am supposed to report the mean value of this signal along with some estimate how accurate this mean value is from a purely statistical point of view. The Idea is: If I conduct the same experiment several times, I get different signals with different mean values. These mean values would be normally distributed. The mean value and standard deviation of this distribution are unknown.

Now if I could perform the experiment several times, estimating the variance of the mean would be straightforward. Then I could estimate some confidence intervals using Student's t-distribution. Yet since compute resources are scarce, this is not an option. All I have is this one sample.

tl;dr:
Is there a way to estimate a confidence interval for the mean when all you have is one sample?

I was thinking about making use of the fact that the signals autocorrelation function vanishes rather quickly. But I haven't been able to turn this into a reliable method.


texts:
I have this extremely distorted and stained image

Is it possible to remove this stain ? Could image inpainting help ?

Please help

EDIT : Another image

After applying anisotropic diffusion and representing image with imagesc(MATLAB)

I tried inpainting however the result isn't good enough

Is there anyway i could improve this output?


texts:
Let's suppose I have finished studying the first college level DSP textbook and I'm looking for a book, that is the next step beyond:
"Digital Signal Processing: Alan V. Oppenheim"

What's the best next book for "Advanced DSP" that would be best to study if on a self-study course?

(self-study meaning, it has sufficient examples and enough technical depth beyond an introduction textbook.)


texts:
I want to do a convolution from the bottom right and not as usual from the top left. I think conv2 of Matlab only does from the top left.

How can I do a convolution in Matlab from the bottom right?

Thank you very much for the answers.


texts:
Can someone give my a link to study to solve this problem.
I can not find any formula for this problem...

Suppose we filter singals with the below pre-emphasis filter:

$y(n)=x(n) - ax(n-1)$

How can I compute $a$ if I know that:

fs=10.000Hz
length 100 in which 98 is zero.
Frequency response at 2100Hz is 1.3429 dB

I don't want the solution only a hint I found this in my notes: 10log10(a^2*w^2+1)=db but with this formula i don't use some information which given.

Thanks


texts: Frequency and amplitude estimation with missing samples

texts:
My input signal is given by $$ x[n] = A\sin(2\pi f n/F_s + \varphi)$$
The signal is uniformly sampled at a sampling rate of $F_s $[Hz]. However some of the samples are missing. Assume that there are not many missing samples and that their location is known (i.e. I know the timestamp for each sample),  I am interested in estimating $A$ and $f$.

My initial approach is to interpolate the vector and use standard methods (e.g. FFT). Another thought is to perform some kind non-linear regression (which I am unfamiliar with).

Any help would be appreciated.


texts: Why is this DFT of a real symmetric signal resulting in complex valued coefficients?

texts:
I am trying to understand exactly how sampling the DTFT to get the DFT works. The signal I'm trying to analyze is $x(n)$ seen below.

$$x(n) = \delta(n\pm2) + 2\delta(n\pm1) + 3\delta(n)$$

Taking the DTFT, we have
\begin{align}
X(\omega)&= \sum\limits_{i=-\infty}^\infty x(n)e^{-j\omega}\\
&= \left(e^{j2ω} + e^{-j2ω}\right) + 2\left(e^{jω} + e^{-jω}\right) + 3\\
&= 2\cos(2ω) + 2(2\cosω) + 3\\
&= 3 + 4\cosω + 2\cos(2ω)
\end{align}

I next implemented this in MATLAB:

x = [1 2 3 2 1];
N = size(x,2);
w = -pi:0.01:pi;
X_DTFT_computational = freqz(x,1,w);
X_DTFT_analytical = 3 + 4*cos(w) + 2*cos(2*w);

This results in the following graphs:

Next, I compute the DFT in two ways:
1. I apply an FFT to the original signal $x(n)$.
2. I sample the DTFT.

The code that achieves this is as follows:

% FFT of x
X_DFT_computational = fftshift(fft(x))

% DFT = Sampled DTFT
X_DFT_analytical(1) = 3 + 4*cos(-4*pi/N) + 2*cos(2*(-4*pi/N));
X_DFT_analytical(2) = 3 + 4*cos(-2*pi/N) + 2*cos(2*(-2*pi/N));
X_DFT_analytical(3) = 3 + 4*cos(0)       + 2*cos(2*(0));
X_DFT_analytical(4) = 3 + 4*cos(2*pi/N)  + 2*cos(2*(2*pi/N));
X_DFT_analytical(5) = 3 + 4*cos(4*pi/N)  + 2*cos(2*(4*pi/N));

Plotting the DFT (while showing the DTFT for comparison), I have

The magnitude of the FFT does give me the correct result.
But just plotting the FFT alone gives me a complex signal (the graph shows the real part).

Why is this signal giving a complex FFT when $x(n)$ is both real and symmetric?


texts: How to achieve full number of words in a Bag-of-Features model based on Hierarchical K-Means?

texts:
I'm implementing a baseline system that uses a Vocabulary Tree, i.e. a BoF model based on HKM, for images classification. Currently I'm obtaining low quality recognition due to the poor quality of the quantization structure resulting from the HKM.

I'm obtaining 100.000 words in my final vocabulary in a tree of depth 6 and branch factor 10 where the theoretical number of words is 10^6 = 1'000.000

In several papers, like those of Zisserman related to Large Scale Landmark recognition, they claim to be using a 1'000.000 words vocabulary, something I found difficult to understand since this number is theoretical while in practice there is no guarantee to obtain it.

Am I understanding something wrong? If not, what should I do to increase the vocabulary size despite of using more descriptors to train the tree?

PS: my only hint so far is using a different seeding algorithm for clusters initialization like k-means++ or gonzales

Thanks, any suggestion is welcome.

After reading @penelope's answer I came back to report my findings. The configuration of my implementation for building the vocabulary tree is: depth 6, branch factor 10, maximum 10 iterations, and random seeding.

Addressing some of the spots:

Severely underpopulated clusters at early levels: I checked clusters at first level and they were kind of uniformly distributed, e.g. the first cluster was typically  10% of the whole set.
K-means termination criteria: I'm running the clustering process until convergence or until reaching 10 iterations. I took this parameter from the experiments of Neister 2006 where the performance will not improve farther than ~12 iterations while doing so could affect the overall algorithm duration.
Uniform coverage of feature space: I am seeding centers using random selection, as I said in the question using other algorithm like k-means++ or gonzalez would become in a more uniform coverage of the feature space.
Number of training features should be considerably bigger than the expected number of clusters: with "expected number of cluster" I understand "expected number of words", being it so, this is very likely the root cause of my problem. I'm training a 1 million words vocabulary using 3 million descriptors while the experiments of Zisserman use ~17 million descriptors from the Oxford Buildings dataset and those of Marc Pollefeys team use ~16 million from the San Francisco Landmark dataset.
Pairing several detectors and descriptors: I have the restriction of comparing my approach which uses binary features against a baseline using DoG/MSER and SIFT descriptors. In the near future I plan to evaluate performance as well using other detectors like AGAST which produce more keypoints of similar quality to those produced by DoG.
HKM termination conditions: in @penelope's answer it was mentioned that obtaining a vocabulary tree of expected size was not a problem. I was wondering if such a thing was truly possible and started to think about HKM termination conditions. In my implementation they are:

Reaching last level.
Having less data than clusters.
Obtaining, at clustering initialization time, centers which are too close one another so there is a big chance they represent the same word, hence not very distinctive.

The first two conditions are enough and necessary but the third one isn't, I included it because the vocabulary might result more distinctive, at least in theory. I borrowed it from the hierarchical k-means index implementation in OpenCV. Could it be convenient removing it?